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ABSTRACT 


In April 2019, the Event Horizon Telescope (EHT) Collaboration reported the first-ever event-horizon-scale images of a black hole, 
resolving the central compact radio source in the giant elliptical galaxy M 87. These images reveal a ring with a southerly brightness 
distribution and a diameter of ~42 uas, consistent with the predicted size and shape of a shadow produced by the gravitationally 
lensed emission around a supermassive black hole. These results were obtained as part of the April 2017 EHT observation campaign, 
using a global very long baseline interferometric radio array operating at a wavelength of 1.3 mm. Here, we present results based on 
the second EHT observing campaign, taking place in April 2018 with an improved array, wider frequency coverage, and increased 
bandwidth. In particular, the additional baselines provided by the Greenland telescope improved the coverage of the array. Multiyear 
EHT observations provide independent snapshots of the horizon-scale emission, allowing us to confirm the persistence, size, and 
shape of the black hole shadow, and constrain the intrinsic structural variability of the accretion flow. We have confirmed the presence 
of an asymmetric ring structure, brighter in the southwest, with a median diameter of 43.3*!2 uas. The diameter of the 2018 ring is 
remarkably consistent with the diameter obtained from the previous 2017 observations. On the other hand, the position angle of the 
brightness asymmetry in 2018 is shifted by about 30? relative to 2017. The perennial persistence of the ring and its diameter robustly 
support the interpretation that the ring is formed by lensed emission surrounding a Kerr black hole with a mass ~6.5 x 10? Mo. The sig- 
nificant change in the ring brightness asymmetry implies a spin axis that is more consistent with the position angle of the large-scale jet. 


Key words. accretion, accretion disks — black hole physics — gravitation — galaxies: active — galaxies: individual: M 87 — 


galaxies: jets 


1. Introduction 


Black holes are a fundamental prediction of Einstein's theory 
of general relativity. A defining feature of black holes is their 
event horizon, beyond which the escape velocity exceeds the 
speed of light, resulting in a dark compact object (Schwarzschild 
1916). Furthermore, photons approaching non-rotating black 
holes within a critical impact parameter Re = V27GM/c?, 
where G is the gravitational constant, M is the black hole 
mass, and c is the speed of light, significantly bend in their 
trajectories by large angles and cannot escape to infinity (e.g., 
Hilbert 1917; von Laue 1921; Bardeen 1973). The value of R, 
changes for spinning black holes, but within <4%. The photon 
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capture and redshift effects produce gravitational signatures in 
the observed images of astrophysical black holes immersed in 
background radiation fields, consisting of a characteristic dark 
center surrounded by ring-like emission (e.g., Luminet 1979), 
commonly referred to as the black hole shadow. Therefore, 
observing the shadow, measuring its properties, such as the 
ring diameter, and comparing it with theoretically predicted 
shadow morphologies provide unique opportunities to explore 
curved spacetime and extreme gravitational potential, allowing 
direct tests of general relativity. This is complementary to tests 
from gravitational waves, high precision astrometry, and cosmol- 
ogy (e.g., Event Horizon Telescope Collaboration 2022f; Abbott 
et al. 2016; GRAVITY Collaboration 2018; Ferreira 2019). 
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From an astrophysical point of view, black holes can be 
formed in many ways, and the masses of astrophysical black 
holes vary greatly, from relatively small (100 Mo; stellar mass) 
to supermassive scales (~10°-10° Mo); for example, readers 
can refer to Greene etal. (2020). Supermassive black holes 
(SMBHs) are predominantly located at the centers of galaxies, 
with larger black holes associated with more massive host galax- 
ies (Kormendy & Ho 2013). When surrounding matter accretes 
onto the central SMBHs, the gravitational potential energy can be 
released in various forms, making active galactic nuclei (AGNs) 
shine across the electromagnetic spectrum, depending on the type 
and rate of mass accretion. Accreting SMBHs can also produce 
collimated relativistic plasma jets and release enormous amounts 
of energy into space in the form of strong magnetic fields and rel- 
ativistic particles. Theoretical studies suggest the ultimate energy 
source of the jet to be either the accretion flow (Blandford & Payne 
1982) or the central black hole itself (Blandford & Znajek 1977), 
in each case mediated through magnetic fields threading the 
plasma in the disk or near the horizon. Therefore, high-resolution 
imaging of accreting SMBHs can also provide unique constraints 
on modeling the detailed physics of the accretion disk and jet 
launching in AGNs (e.g., Blandford et al. 2019). 

Since the angular diameter of the black hole shadow is 
dg = 2R.D! =~ 10GMc~?D7! where D is the distance 
to the black hole (Bardeen 1973), the imaging capability of 
the black hole shadow requires a mass-to-distance ratio suffi- 
ciently large to provide dsn resolvable with the existing instru- 
ments. In this regard, SMBHs in the center of the radio 
galaxy M 87 (hereafter M 87*) and our Galaxy, Sagittarius A* 
(Sgr A*), offer the best opportunities. As for M87, its dis- 
tance is measured via multiple methods using primary or sec- 
ondary distance indicators, and we adopt D = 16.8 + 0.8 Mpc 
(Bird et al. 2010; Blakeslee et al. 2009; Cantello et al. 2018, 
see Event Horizon Telescope Collaboration 2019f). The central 
black hole mass is approximately 6 x 10? Mọ (Gebhardt et al. 
2011; Event Horizon Telescope Collaboration 2019f). The com- 
bination of D and M yields a d,, ~ 40 pas, which is the second 
largest value on the sky after Sgr A* (M « 4 x 106 Mo, D « 
8 kpc, da, & 50 was; see Event Horizon Telescope Collaboration 
2022a) M87 is a representative Fanaroff-Riley I galaxy 
(Fanaroff & Riley 1974) with a prominent jet extending 
from the vicinity of the central black hole to kiloparsec 
scales (Owen et al. 1989; Reid et al. 1989; Junoret al. 1999; 
Asada & Nakamura 2012; de Gasperin et al. 2012; Hada et al. 
2016; Mertens et al. 2016; Kim etal. 2018a; Nakamura et al. 
2018; Walker et al. 2018; Park et al. 2019; Goddi et al. 2021; 
Lu et al. 2023). Multi-wavelength (MWL) studies of the jet from 
radio to y-rays show both steady and sporadic bright flaring 
emission (see EHT MWL Science Working Group 2021). This 
makes M 87 one of the best candidates for ultra-high resolution 
imaging of an astrophysical black hole and the jet launching site. 

The Event Horizon Telescope (EHT) is a global net- 
work of radio telescopes for very long baseline inter- 
ferometry (VLBI) observations primarily at a wavelength 
of A x 1.3mm or a frequency of v « 230GHz 
(Event Horizon Telescope Collaboration 2019b). In VLBI, we 
measure complex visibilities or Fourier components of the 
radio brightness distribution of the sky at spatial frequencies 
u = (u,v), corresponding to the projected baselines in units 
of observing wavelength in a plane normal to the direction 
of the phase reference position (Thompson et al. 2017). These 
complex visibilities can be separated into amplitude and phase 
components, and are commonly constructed as closure quanti- 
ties (Rogers et al. 1974; Readhead et al. 1980; Blackburn et al. 


2020), which are immune to station based systematic corrup- 
tions. The longest baselines of the EHT result in an array with 
a theoretical diffraction limited resolution of 1/ |u| ~ 25 uas 
from the ground, and can directly resolve the black hole shad- 
ows in M87* and Sgr A*. The first image of the black hole 
shadow of M 87* was obtained by the EHT observations in April 
2017 (Event Horizon Telescope Collaboration 2019a,b,c,d,e,f, 
2021a,b, 2023, hereafter M 87* 2017 I-IX). Also, another image 
of the black hole shadow, that of SgrA*, was published 
(Event Horizon Telescope Collaboration 2022a,b,c,d,e,f, here- 
after Sgr A* 2017 I- VI). 


2. Previous EHT results and outline for new results 


Below we provide a brief review of the major results on M 87* 
from the 2017 EHT observations to introduce the 2018 EHT 
results, which are the main focus of this paper. For additional 
details, readers can refer to M 87* 2017 I and Sgr A* 20171 for 
an overview of the first black hole shadow imaging of M 87* and 
Sgr A*, respectively. 


2.1. Results from previous EHT observations of M 87* 


In April 2017, eight radio telescopes executed the first EHT 
VLBI observations of M87* at a wavelength of 1.3mm 
(M 87* 2017 ID). Three independent EHT data calibration 
pipelines were developed and applied to validate the fringe 
detection and quantify their systematic uncertainties to be used 
for downstream analysis (M 87* 2017 III). The M 87* data reveal 
the presence of two nulls in correlated flux density at ~3.4 
and ~8.3 GA (M 87* 2017 III) and temporal evolution in closure 
quantities, indicating intrinsic variability of compact structure 
on a timescale of days or several light-crossing times for a few 
billion solar-mass black hole (M 87* 2017 IV). 

Independent analysis of the calibrated data using imaging 
and geometric modeling methods produce reconstructions that 
are well characterized by a bright and asymmetric ring-like 
structure with a diameter of 42 + 3 uas, a central flux depres- 
sion with an intensity ratio of >10:1, and a position angle of the 
brightness asymetry of ~180°. These measurements are consis- 
tent with expectations for a SMBH of mass M = (6.5 + 0.7) x 
10? Mo (M 87* 2017 IV; M 87* 2017 VI). A small level of inter- 
day variability was also detected, which was analyzed in follow- 
up studies (e.g., Georgiev et al. 2022; Satapathy et al. 2022). 
Analyses from outside the EHT collaboration featuring both 
imaging (Carilli & Thyagarajan 2022; Arras etal. 2022) and 
geometric modeling methods (Lockhart & Gralla 2022) agree 
that the 2017 EHT data are consistent with a ring-like structure. 
There is another external analysis of the 2017 data that claims 
not to support the EHT interpretation (Miyoshi et al. 2022). 

In addition, Wielgus et al. (2020) reported results of 1.3 mm 
VLBI monitoring of the compact core of M 87 with a sparser 
pre-EHT array. By assuming a ring-like structure, this analysis 
found a consistent diameter across observations from 2009 to 
2017 and strong evidence that the position angle of the bright- 
ness asymmetry varies by many degrees from year to year. 
Linear polarized features were also detected with a ring-like 
structure (M 87* 2017 VID, along with the first limits on the cir- 
cular polarization at the horizon scales (M 87* 2017 IX). 


2.2. Summary of the new results and paper organization 


The EHT array in 2018 was significantly improved since 2017. 
The recording rate has increased by a factor of two compared 
to 2017, and the Greenland Telescope (GLT; Inoue et al. 2014; 
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Fig. 1. Representative example images of M 87* from the EHT observations taken on 2017 April 11 and 2018 April 21 (north is up and east is to the 
left). The 2017 image is generated with the average of fiducial parameter sets from the imaging techniques used in M 87* 2017 IV. The 2018 image 
is created by taking the average of the blurred images generated by the imaging techniques found in Sect. 5. Comparison of the images shows 
consistency in the diameter across observation epochs, but a shift in position angle of brightness asymmetry. The circle represents a Gaussian 


blurring kernel with a full width half maximum of 20 uas. 


Chen et al. 2023) joined the EHT array for the first time. The 
addition of this new station results in a better (u,v) coverage 
compared with the 2017 EHT observations, particularly along 
the north-south direction. New EHT images of M 87* can pro- 
vide further evidence that the event-horizon-scale images from 
the 2017 observations are consistent with the prediction from 
general relativity that strongly lensed emission around a black 
hole should form a persistent ring-like structure. A detection of 
variability around the ring-like structure can also provide us with 
a better insight into the dynamics of the underlying plasma and 
the black hole angular momentum. 

The structure of this paper and the main results are as fol- 
lows: In Sect. 3, we describe the observations and processing of 
the 2018 EHT campaign data, providing a summary of the data 
properties and issues compared to that from the 2017 EHT cam- 
paign. As in 2017, we see a clear indication of a deep null in 
the visibility amplitudes around 3.4 GA. We also see significant 
closure phase deviations relative to 2017, indicating non-trivial 
structural changes. 

We explain how we derived pre-imaging constraints on the 
expected compact flux density and its size in Sect. 4. Combin- 
ing results from simple visibility modeling and information from 
lower frequencies, we settled on pre-imaging initial parameters 
for the compact flux density to be between 0.4 to 1.0Jy and the 
size of the source to be <100 pas. 

In Sect. 5, we describe the imaging procedure and investi- 
gate the primary image morphologies. We reconstruct the images 
using a variety of imaging techniques, ranging from traditional 
inverse imaging to Bayesian posterior exploration. We compare 
the recovered image structure across the different imaging meth- 


A79, page 4 of 63 


ods, frequency bands, calibration pipelines, and observing dates. 
We recover a clear ring-like structure, with a deep central depres- 
sion and a diameter of ~42 uas. We also find that the position 
angle of the brightest region of the ring is between 200° and 
230°, which is different from the position angle measured in 
2017. We present the average images of M 87* across methods 
on 2018 April 21 compared to 2017 April 11 in Fig. 1. 

Inspired by the clear ring-like structure seen in the imag- 
ing methods, in Sect. 6 we quantify the support for the ring-like 
structure by comparing the Bayesian evidence of different geo- 
metric models against the data. We find that ring-like models are 
strongly preferred by the 2018 data. We then check the consis- 
tency of important physical quantities such as the diameter and 
position angle of the emission ring by directly modeling those 
features in the visibility domain. We find that all direct modeling 
methods find a position angle around 210°. The direct modeling 
methods also prefer a slightly higher diameter, around 45 pas. 
The diameter estimates from imaging and direct modeling live 
within each other's error bars, and systematic differences in the 
median diameters can be attributed to model-dependent resolu- 
tion effects. 

In Sect. 7, we build upon the findings from the 2017 EHT 
campaign by comparing them against the imaging and direct 
modeling results from the 2018 data. We explore the persis- 
tence of the ring diameter between 2017 and 2018 and discuss 
the stability of the mass across the two years. We also discuss 
the possible origin of the shift in the position angle between the 
two years, which is consistent with our understanding of turbu- 
lent material around the black hole. We investigate the compact 
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Fig. 2. Map showing the stations that participated in the EHT 2018 cam- 
paign (black circles), which differs from the EHT array in 2017 by the 
addition of the GLT. Co-located sites in Chile and Hawai ‘i appear super- 
imposed. The SPT projected location from the back of the map is indi- 
cated with a dashed circle, and baselines to this station are represented 
with dashed-lines. While the SPT cannot observe M 87^, it observed 
3C 279 and was used to calibrate the data. 


flux density estimates produced by the Bayesian image recon- 
struction and modeling methods, and find that the compact flux 
density is especially difficult to constrain in the 2018 data. The 
main conclusions of our work are summarized in Sect. 8. A 
detailed theoretical interpretation will be presented in a later 


paper. 


3. Observations and data processing 


In this section, we briefly describe the 2018 EHT observing cam- 
paign (Sect. 3.1), the data correlation (Sect. 3.2), the calibration 
and data reduction (Sect. 3.3), and finally provide some compar- 
isons between the 2017 and 2018 data properties (Sect. 3.4). 


3.1. Overview of the 2018 observing campaign 


The 2018 EHT observing campaign was scheduled from 2018 
April 18 to 29. Based on the expected opacity and weather con- 
ditions at each site, observations were triggered on six observ- 
ing days with a total of nine participating stations, including 
all the stations which participated in the 2017 EHT campaign 
(see M 87* 2017 II) with the addition of the GLT. Thus, the 
2018 array includes the phased ALMA and the phased Sub- 
millimeter Array (SMA), made up of 43 12-m and seven 6- 
m dishes, respectively, as well as the following seven single- 
dish stations: Atacama Pathfinder Experiment (APEX), IRAM 
30-m telescope (PV), South Pole Telescope (SPT), Large Mil- 
limeter Telescope Alfonso Serrano (LMT), Submillimeter Tele- 
scope (SMT), James Clerk Maxwell Telescope (JCMT), and the 
GLT. In Fig. 2, we show a map of the EHT array with the stations 
that participated in the 2018 campaign. 


Table 1. Median zenith sky opacities (1.3 mm) at EHT sites during the 
2018 April observations toward M 87*. 


Station Median zenith T1.3 mm 
Apr.21] Apr.22  Apr25  Apr.28 

ALMA 0.10 0.10 0.12 - 
APEX 0.10 0.12 0.14 - 
GLT 0.25 0.41 0.43 0.29 
LMT 0.08 0.20 - = 
SMT 0.21 0.16 0.32 0.50 
JCMT 0.24 - 0.09 0.36 
PV 0.39 - 0.44 0.18 
SMA 0.17 - 0.07 0.35 


Four frequency bands centered at sky frequencies of 213.1 
(band 1), 215.1 (band 2), 227.1 (band 3), and 229.1 GHz (band 4) 
were recorded, except for the GLT which observed only in bands 
3 and 4. Each band has a bandwidth of 2048 MHz for every sta- 
tion, except for ALMA, which provides an effective bandwidth 
of 1875 MHz. This represents an improvement in spectral cover- 
age by a factor of two compared to the 2017 observations. Addi- 
tionally, the LMT had an increased sensitivity, given its increase 
in effective dish diameter from 32.5 m in 2017 to 50m in 2018 
and a change from a double-sideband (DSB) receiver to a two- 


single-sideband (2SB) one. 
The telescope systems recorded both right-hand and left- 


hand circular polarization (RCP and LCP), except for JCMT, 
which observed only RCP throughout the 2018 campaign, and 
ALMA, which recorded linear polarization. Details on the treat- 
ment of different polarization types are given in Sect. 3.2. 

Scans on M 87* were included in four observing days with up 
to eight participating stations. The scan duration varied between 
4 min and 7 min through the observations (although shorter scans 
were used in two days of the campaign). For calibration pur- 
poses, M 87* scans were interleaved with 3—4 min duration scans 
on the source 3C 279. The median zenith sky opacities at 1.3 mm 
throughout the array on each observing day are provided in 
Table 1. In Fig. 3, we show the scan distribution per participating 
station in the schedule of the best observing days toward M 87* 
and 3C 279. 

The observations on April 28 had an overall higher median 
zenith opacity throughout the array than the other observing days 
(Table 1) and were notably affected by bad weather at ALMA. 
This prevented the detection of fringes in most baselines and 
resulted in many triangles with no closure information. Conse- 
quently, observations on April 28 did not pass the required qual- 
ity checks (including the second phase of the ALMA quality 
assurance process, i.e., QA2) and were discarded in subsequent 
analyses. 

The operations at the LMT were halted midway through the 
2018 EHT campaign owing to local security issues. On the first 
day of the campaign (April 21), the observations started later 
than planned due to technical issues, which repeated in the mid- 
dle of the observing day. At the start of that same day, the first 
few scans from Hawaiʻi stations were lost to bad weather, and 
the SMA started observations even later due to technical issues. 
This resulted in less observing time for the LMT and Hawaiʻi 
stations, which explains the differences in the (u, v) coverage in 
2018 compared to 2017 (Fig. 4, left panels). 

Furthermore, the (u,v) coverage on April 22 is minimal. 
The observations on April 25 suffered from coherence loss on 
baselines with APEX, bad weather at PV, and poor phasing 
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Fig. 3. EHT observing schedules for M 87* (blue) and its calibrator 3C 279 (orange) on the 2018 April 21 (left panel) and April 25 (right panel) 
observing days, which began at the end in UTC of April 20 and 24, respectively. Open rectangles represent scans that were scheduled but not 
observed owing to weather or technical issues. The filled rectangles mark the scans with detections in the final data set. On these two particular 
days, scan durations are typically 4 to 5 min each for M 87* and 4 min each for 3C 279, as reflected by the width of each rectangle. 


efficiencies at the SMA. Hence, April 21 is the best observing 
day of the 2018 EHT campaign for M 87*, followed by April 25. 


3.2. Correlation 


The data from the 2018 observing campaign were correlated at 
two correlation centers located at the Max-Planck-Institut für 
Radioastronomie in Bonn, Germany and at the MIT Haystack 
Observatory in Westford, Massachusetts, USA. The Distributed 
FX (DiFX, version 2.6.2) correlation package (Deller et al. 
2011) was used. The data underwent multiple correlation passes 
to diagnose and correct data issues, including incorrect polar- 
ization labeling on the GLT and LMT. The data used in this 
paper is based on the fourth revision (Rev4) of the 2018 cor- 
relation products, released as part of the 2018 EHT Level 1 (L1) 
data package that also includes metadata required for absolute 
flux density calibration (Koay et al. 2023a). The final correlation 
produced 32 baseband channels, each 58 MHz wide with a spec- 
tral resolution of 0.5 MHz and averaged to a 0.4 s accumulation 
period. After correlation, the recorded RCP and LCP streams 
are multiplied to create the parallel-hand (RR, LL) and cross- 
hand (RL, LR) polarization products. Due to the non-matching 
sampling rates between ALMA and the other stations, a por- 
tion of the bandwidths is lost during correlation, resulting in 
a correlated bandwidth of 1.856 GHz per band (Matthews et al. 
2018, M 87* 2017 III). All the EHT stations observe in circular 
polarizations except ALMA, which observes in linear polariza- 
tion. The PolConvert v1.7.9 (Martí-Vidal et al. 2016) was then 
run to convert the mixed polarization data products (XL, XR, 
YL, YR) to circular polarization basis on baselines to ALMA, 
described in more detail by Goddi et al. (2019). 


3.3. Calibration and post-processing 


To maximize fringe amplitudes when averaging the data in 
time and frequency, residual delays, delay rates, and atmo- 
spheric phase turbulence need to be corrected. Additionally, 
sampling losses and telescope bandpasses need to be calibrated. 
These signal stabilization (Janssen et al. 2022) data reduction 
steps are performed with two independent pipelines, EHT-HOPS 
(Blackburn et al. 2019) and the CASA-based (CASA Team 
2022; van Bemmel et al. 2022) rPICARD software (Janssen et al. 
2019a). These pipelines are described in more detail in 
M87* 201710 and Sgr A* 2017II. No significant algorithmic 
changes were made to the EHT-HOPS pipeline fringe fitting 
methods compared to the 2017 analysis. For the rPICARD 
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pipeline, we now use a priori determined solutions for the typ- 
ically ~7.3 ns single-band and multi-band delay offsets of the 
phased ALMA instead of solving for these offsets using the 
ALMA-APEX baseline as was done for the 2017 data. 

The stabilized data are then loaded into the AIPS software 
package (Greisen et al. 2003), where a priori amplitude cali- 
bration is applied using the task APCAL to convert the corre- 
lation coefficients to a common flux density scale across the 
very heterogeneous EHT array. These are based on calibra- 
tion tables in the standard ANTAB format containing informa- 
tion on telescope sensitivity, elevation-dependent gain curves, 
and time-dependent system temperatures derived from the raw 
metadata collected for each station, as described in Koay et al. 
(2023a). A summary of the amplitude calibration parameters 
of each station and their associated uncertainties is provided in 
Appendix A. The data are then averaged over time (10s) and 
each 1.856 GHz band before being exported by AIPS as UVFITS 
files. Subsequently, we perform a polarimetric gains ratio cali- 
bration (polgainscal) to align the RR and LL phases in the 
HOPS data. We correct for phase offsets between the two polar- 
ization streams and allow the Stokes 7 visibilities to be obtained! . 
Network calibration (M 87* 2017 IIT) with 10 s solution intervals 
is then applied to correct the amplitudes of stations with a co- 
located partner (ALMA and APEX in Chile, SMA and JCMT 
in Hawai‘i), using redundancies and total flux density measure- 
ments from the phased ALMA (EHT MWL Science Working 
Group et al., in prep.). The network calibration takes advantage 
of the fact that the source is dominated by an unresolved core 
component on the intra-site baselines and that the flux densi- 
ties on the intra-site baselines are comparable to that of the total 
flux densities of the core component as observed by the phased 
ALMA. We expect that the difference between the correlated 
flux density measured by the phased ALMA and the ALMA- 
APEX baselines should be less than 1%. 

At each stage of the calibration, the data products are veri- 
fied via consistency checks of closure quantities between bands 
and polarizations, as well as checks of the trivial closure quan- 
tities (Fig. B.1, see also M 87* 2017 III and Wielgus et al. 2019, 
for more details). We also use these closure analyses, and the 
more novel closure trace analyses (Broderick & Pesce 2020) to 
roughly quantify the non-closing systematic uncertainties pos- 
sibly arising from polarimetric leakages and bandpass effects 


! Tn the CASA data, a priori corrections for the field rotation angle vari- 
ations are applied upstream and the RCP-LCP phase offsets are solved 
as part of the instrumental calibration steps. 
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Fig. 4. M 87* (u,v) coverage (colored points) in band 3 (top panels) and band 2 (bottom panels) for observations on 2018 April 21 (left panels) 
and April 25 (right panels), overlaid on the low-band (u,v) coverage for 2017 April 11 (gray points). The dashed circles show baseline lengths 
corresponding to fringe spacings of 25 and 50 uas respectively. The (u, v) coverage in band 1 and band 4 is comparable to that in band 2 and band 


3, respectively. 


(see Appendix B). For M 87*, the non-closing systematic errors 
are estimated to be 1.2? in closure phases and 2.7% in log clo- 
sure amplitudes, using both the HOPS and CASA data products. 
Assuming the errors are baseline independent, these translate to 
0.7? systematic errors in visibility phases and 1.3% systematic 
error in visibility amplitudes. These uncertainties do not include 
station-based amplitude gain errors such as those described in 
Appendices A and E, which are corrected by model-dependent 
calibration during imaging. 

Since some imaging pipelines use time-averaged data up to 
the length of a single scan, that is 4 to 5 min, we examine the level 


of decoherence losses in amplitude when the data are averaged in 
time. The details are presented in Appendix C. Such amplitude 
losses arise from imperfectly aligned phases after calibration and 
can lead to non-closing errors in visibilities. We find that 92.9% of 
HOPS and 85.6% of CASA M 87* data have >90% coherence (i.e., 
exhibit decoherence losses of less than 1096) when averaged over 
the entire scan. When considering April 21 data only, the coher- 
ence levels are better, with 97.4% and 87.9% of HOPS and CASA 
scan-averaged data having better than 9096 coherence. 

We also check for consistency between the closure phases 
and log closure amplitudes derived from the HOPS and CASA 
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Fig. 5. Measured correlated flux densities of M 87* as a function of baseline lengths in units of wavelength, for 2018 April 21 (top panels) and 
April 25 (bottom panels) in band 3, for both HOPS (left panels) and CASA (right panels) outputs. The 2018 data (colored points) are overlaid 


on the corresponding flux densities of the 2017 April 11 observations 


in low-band (gray points). All data shown include a priori station-based 


amplitude calibration and network calibration but are prior to any model-dependent self-calibration. Error bars denote +10 from thermal noise. 
Redundant baselines are shown with different symbols: circles for baselines to ALMA and SMA; diamonds for baselines to APEX and JCMT. 


The dashed line corresponds to an azimuthally symmetric thin ring mod 


data products for which signal-to-noise ratio (S/N) >7 and find 
that there are <4% outliers of >20 when o is equivalent to 
the thermal noise (om). The fraction of outliers decreases to 
<2% when o includes the systematic uncertainties described 
above. In Fig. B.1, the fraction of 23c' outliers are consistent 
between HOPS and CASA also indicate that both pipelines pro- 
duce consistent calibrated data sets. For reference, we provide 
detailed discussions of data issues in Appendix D, and in partic- 
ular, highlight gain calibration issues for the GLT and LMT in 
Appendix E. 


3.4. M87* 2018 data properties and comparisons with 2017 


The 2018 EHT observations were made with an array similar to 
the one in 2017, except for the addition of the GLT and dou- 
bling the bandwidth of the 2017 observations, as mentioned in 
Sect. 3.1. The central frequency of bands 3 and 4 in 2018 coin- 
cides with the low- and high-bands in 2017, respectively. In 
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el with a 42 uas diameter. 


Fig. 4, we show the (u, v) coverage of the best observing days 
in 2018 (April 21 and 25) in bands 2 and 3, overlaid on the (u, v) 
coverage of the best observing day in 2017 (April 11) in the low- 
band. 

The GLT adds north-south coverage to the EHT array in 
bands 3 and 4, probing baseline lengths of ~7.1 GA with the 
co-located Chilean sites and baseline lengths of ~3.8 GA and 
~2.8 GA with intermediate stations. However, on April 21, the 
longest baselines in the east-west direction given by PV-Hawai'i 
are missing, such that the east-west coverage is worse than that 
of April 25 and the 2017 observations. The LMT is missing on 
April 25, and its coverage is significantly reduced on April 21 
compared to 2017, which impacts the compact flux density con- 
straints in 2018. 

Figure 5 shows the visibility amplitude versus (u, v) distance 
for 2018 April 21 and 25, overlaid on 2017 April 11. The a priori- 
and network-calibrated visibility amplitudes (see M 87* 2017 HI, 
for details) display the characteristic secondary peak beyond 
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Fig. 6. Variation of closure phases as a function of Greenwich mean sidereal time (GMST) for selected triangles using HOPS calibrated data. Red 
diamonds denote data from 2017 April 6 (low-band), gray circles denote data from 2017 April 11 (low-band), and blue squares denote data from 


2018 April 21 (band 3). Error bars are the +1o uncertainties. 


a deep amplitude minima at —3.4 GA, also observed in 2017 
(M 87* 2017 IV), despite the differences in the (u, v) coverage. 
A ring-like structure will also produce a second null beyond 
~8.3 GA. However, the longest baseline the EHT can probe is 
8 GA between the Hawai‘i and PV stations, which has only one 
detection on 2018 April 25. While the S/N at the amplitude min- 
ima is lower (x10) than that at the other regions, they provide 
very strong upper limits on the amplitudes (x50 mJy) at these 
"null" locations. In particular, part of the LMT-Hawai'i cover- 
age samples similar fringe spacings as the GLT-PV and GLT- 
SMT baselines (although at almost perpendicular orientations), 
and part of the Chile-LMT coverage samples similar spacings 
as the GLT-LMT baseline. The GLT-PV baselines display lower 
amplitude than the GLT-LMT and GLT-SMT baselines. Thus, in 
2018 we also see evidence of source anisotropy around the first 
null thanks to the GLT, as was observed in 2017 when the GLT 
was not available, but the corresponding fringe spacings were 
sampled by the 2017 LMT baselines. 

The closure phases in all non-trivial triangles show dif- 
ferences in 2018 compared to 2017, providing evidence for 
significant changes in structure between these two epochs. In 
Fig. 6, we compare 2017 and 2018 coherently averaged visibil- 
ities on the same triangles (or similar ones) selected in 2017, 
where day-to-day variability of closure phases was observed 
(see Fig. 14 in M 87* 2017 III). Here we also include a wide 
GLT triangle whose closure phases deviate from zero, indi- 
cating the presence of resolved asymmetric structure, similar 
to what has been observed in other triangles in 2017 before 
the inclusion of the GLT in the EHT array (M 87* 2017 III). 
Gray circles and red diamonds in Fig. 6 show closure phase 
measurements from two days in 2017, and blue squares show 
closure phases measured in 2018. While some triangles show 
some small deviations in closure phase between days in 2017, 
the 2018 closure phases are qualitatively different, demonstrat- 
ing clear indications of different asymmetric structure in the 
2018 data. 


4. Compact flux density and source size constraints 


In this section, we analyze the network-calibrated visibilities 
of M87* to produce estimates of the compact flux density and 
source size. These are used to generate reasonable synthetic data 
sets (Sect. 5) and inform our choices of imaging parameters 
over which to survey for the non-Bayesian image reconstruction 
methods (see Sect. 5). 

Measurements of the total flux density of the arcsecond- 
scale radio core in M 87 at 1.3 mm are between 1.0 and 2.0Jy 
(Bower et al. 2015). Previous VLBI observations of M 87 at 
1.3 mm between 2009 and 2017 indicated total compact flux 
densities at the milliarcsecond and microarcsecond scale fluc- 
tuating within a range of 0.5-1.0Jy (Doeleman et al. 2012; 
Akiyama et al. 2015; Wielgus et al. 2020). We infer that the dis- 
crepancy between the total emission at arcsecond scales and 
VLBI scales is associated with the unconstrained extended emis- 
sion outside of the compact emission region, presumably domi- 
nated by the jet. 

Angular scales in M 87* are sampled over nearly five orders 
of magnitude in (u,v) space by the EHT baselines. There is 
a large gap in the coverage between the intra-site baselines 
(ALMA-APEX and SMA-JCMT) of 0.1—1.1 MA and the inter- 
site baselines of 1.3—7.3 GA, which are sensitive to the sub- 
arcsecond and microarcsecond scale structures, respectively. The 
longest intra-site baseline is ALMA-APEX, with fringe spac- 
ings of 106—131 mas for M 87*. On the other hand, that of the 
shortest inter-site baseline (SMT-LMT) is 139—166 uas. There- 
fore, the structures on spatial scales of ~0.2—100 mas are not 
sampled in the visibility domain, and are thus unconstrained by 
the current EHT measurements. 

The EHT data provide two relevant estimates of the total 
flux densities at different spatial scales. The first one is the total 
flux density of the arcsecond-scale radio core, constrained by the 
short intra-site baselines. The total flux density as measured by 
the ALMA-only interferometric observations was Fi ~ 1.13 Jy 
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by averaging the values over four bands on April 21 and 22. The 
second one is the total compact flux density, F;y < Frot, which 
is the integrated flux density of the microarcsecond-scale struc- 
ture constrained by inter-site baselines. 

We estimate the range of F;y. and the size of the compact 
emission region pct, represented by an equivalent Gaussian full- 
width half maximum (FWHM), following the two independent 
procedures in M 87* 2017 IV. One is directly estimated from the 
2018 EHT visibility data. With this procedure, we can derive 
the constraints on Fepcet and @cpct simultaneously. The other is 
a constraint on the Foyt by utilizing quasi-simultaneous VLBI 
data based on an extrapolation between 1.3 mm and the longer 
wavelengths. In Appendix F, we derive a series of constraints on 
Fepct and 6,54 using these procedures. 

The EHT constraints give Fp = 0.30-1.13 Jy with a Oper 
between 39 and 98 uas. The source size constraints are roughly 
consistent with the 2017 constraints. The resulting constraints 
on Fepet for 2018 are looser than those for 2017 due to the lack 
of additional constraints using SMT-PV and LMT-PV cross- 
ing baseline tracks. Nevertheless, the MWL constraints give a 
slightly tighter range of Fep = 0.5-0.7Jy within the cen- 
tral 100 x 100 uas? of the nucleus, which is comparable to the 
Fa = 0.56-1.03 Jy estimated from the 2017 data. Therefore, 
we conclude that the total compact flux density of M 87* in 
2018 at 1.3mm is broadly consistent with that of 2017 when 
considering all these estimates. While these flux density and 
size constraints in the pre-imaging analysis provide a valuable 
guide for subsequent imaging, analysis, and synthetic data gen- 
eration, we do not strictly enforce them and leave the com- 
pact flux density as a free parameter in the subsequent imaging 
surveys. 

In our pre-imaging constraints obtained above we pre- 
dict a compact source size less than 100 uas with an uncon- 
strained extended structure larger than ~0.2 mas. Since there 
are no visibilities probing structure above —0.2 mas, the EHT 
data cannot be used to constrain jet emission above this scale. 
Moreover, the emission from the extended jet is significantly 
lower in surface brightness than the core region, even in 
3 mm observations (Hada et al. 2016; Kim et al. 2018a; Lu et al. 
2023), so the jet emission is expected to be even weaker 
at 1.3 mm due to its optically thin (steep) synchrotron spec- 
tral nature (EHT MWL Science Working Group 2021). Even 
though we have improved coverage with the GLT in 2018, 
the continued lack of short baselines means we do not 
expect the 2018 EHT data can constrain this very dim jet 
emission. 


5. Imaging and image domain analysis 


The EHT's sparse (u, v) coverage results in an ill-posed inverse 
problem that prevents the recovery of a unique image from a 
measured set of visibilities. By using realistic priors, it is nev- 
ertheless possible to reconstruct brightness distributions that 
are consistent with the data. While these reconstructed dis- 
tributions are still non-unique, they should represent a fam- 
ily of images that are interpretable within our understanding 
of the physical system. While it is possible to produce rea- 
sonable image reconstructions using only closure quantities, 
all methods featured in this section utilize some data prod- 
ucts that are complicated by persistent systematic calibration 
errors in both the data amplitude and phases. This requires 
each reconstruction method to undergo a careful self-calibration 
or station gain reconstruction process in order to extract the 
most information from the data. As in M 87* 2017IV, both 
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inverse and forward modeling were used in the image recon- 
structions. For inverse modeling, a CLEAN-based algorithm 
(e.g., Hógbom 1974; Clark 1980) was employed with DIFMAP 
software (Shepherd 1997, 2011). The forward modeling tech- 
niques used in this study were regularized maximum likeli- 
hood (RML) methods represented by eht-imaging (Chael et al. 
2018, 2019a,b) and SMILI (Akiyama et al. 2017a,b), and two 
Bayesian sampling methods to explore the full posterior space 
with THEMIS (Broderick et al. 2020a,b, 2022a) and Comrade 
(Tiede 2022). 

To explore how the reconstructed images were affected by 
various imaging and optimization choices, we generated syn- 
thetic data sets using ten geometric models (cres-96, cres, 
cres96, cres180, dblsrc, disk, ring, edisk, point+disk, 
point+edisk) and a snapshot of a general-relativistic magne- 
tohydrodynamic (GRMHD) simulation (see Appendix G). The 
synthetic data were designed to have properties similar to the 
EHT M 87* visibility amplitudes, including prominent ampli- 
tude nulls, but also reflected a diversity of both ring-like and 
non-ring source structures. This procedure was similar to previ- 
ous EHT imaging analyses (M 87* 2017 IV; Sgr A* 2017 IIT). For 
the RML and CLEAN imaging methods, four of these ten data 
sets (cres186, ring, dblsrc, disk) were used as training sets 
to evaluate the ability of unique combinations of imaging param- 
eters to faithfully reconstruct images close to ground truth. For 
the Bayesian methods, all eleven data sets were used as valida- 
tion exercises to confirm the ability of the Bayesian methods to 
reconstruct diverse image structures. 

Section 5.1 presents the image strategies used by each 
method. Section 5.2 presents the images and discusses consis- 
tency between different methods, frequency bands, calibration 
pipelines, and observation dates. Section 5.3 describes the pro- 
cess by which we measure important image domain quantities 
such as the diameter and position angle. 


5.1. Strategy of the imaging analysis 


We conducted parameter surveys with three imaging pipelines; 
employing CLEAN using the DIFMAP software, and RML meth- 
ods using eht-imaging and SMILI. We used the four training 
data sets to explore the impact of different imaging assumptions, 
such as hyperparameters (weights for both the data and regu- 
larizers, see M 87* 2017 IV) and optimization choices, on the 
resulting image morphology. All images in the surveys were 
reconstructed from data calibrated as described in Sect. 3.3. 
From these surveys, we selected a “Top Set" of parameter com- 
binations for each method, which represent the set of best-fit 
images to the data. Each synthetic data set was blurred with a 
nominal beam width corresponding to the effective resolution 
of the longest baseline. We then compared the normalized cross 
correlation, Oxx, of the unblurred ground-truth image against 
the blurred ground truth to determine a pyx cutoff value. We 
then found the model pyx value for each set of image param- 
eters against the blurred synthetic data, and a set of model 
image parameters passed to the next step if the model pyx is 
above the cutoff value described above. From this pruned set, 
the Top Set images were then selected by calculating the nor- 
malized closure quantity chi-squares of the images for the real 
M 87* data, and taking only the images that have a normalized 
X) < 2 on data with 0% systematic noise (to account for non- 
closing errors) for the RML methods and 1046 systematic noise 
for CLEAN (see Sect. 6.3.1 in M 87* 2017 IV and Sect. 5.2.1 in 
this paper for more details). The closure quantities were aver- 
aged over the entire scan before comparison to ensure sufficient 
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Table 2. Parameter survey results for April 21 band 3 data. 


DIFMAP (1260 Param. combinations; 303 in Top Set) 


Compact 0.4 0.5 0.6 0.7 0.8 

flux (Jy) [13% | 17% 24% 239 239 

ALMA 0.01 0.1 0.3 0.5 0.7 1.0 

weight factor 20% [24%] 26% 17% 10% 3% 

Mask diam. 40 50 60 70 80 90 100 
(pas) 0% 4%  [27%| 27% 23% 12% 796 
(u, v) weight -2 -1 0 

exponent K 32% [59% | 9% 

Stop Flux reached Arms x 10* 
condition [52% | 48% 


eht-imaging (12288 Param. combinations; 874 in Top Set) 


Compact 0.4 0.5 0.6 0.7 
flux (Jy) 18% [27%] 29% 269 
Init./MEM 40 50 60 
FWHM (uas) |57%| 39% 4% 
Systematic 0% 1% 2% 5% 
error 14% 22% 30% 34% 
Regularizer: 0 1 10 100 
MEM 0% 0% 5% 95% 
TV 33% |36% 30% 0% 
TSV 31% |33% 35% 1% 
& 2690 |25%] 25% 24% 
SMILI (36288 Param. combinations; 5333 in Top Set) 

Compact 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
flux (Jy) 10% 18% [18%] 18% 15% 13% 10% 
EW Soft mask 40 50 60 70 80 90 100 110 
FWHM (uas) 13% 15% 15% [ 13% | 12% 11% 11% 11% 
Systematic 0% 1% 2% 
error 36% 33% 30% 
Regularizer: 0 10 10° 10° 104 10° 
TV 19% 19% 20% 23% 19% 0% 
TSV 19% 19% 20% 23% 18% 0% 

0 107? 107! 1 10 10° 
w O% 0% 4% 35% |58%| 2% 


Notes. Below each parameter values we specify the fraction of the Top Set parameter combinations that include that value. Boxed parameters are 
those corresponding to the fiducial images. The fiducial parameters are determined by identifying the parameter combinations that jointly perform 
best on all the synthetic data sets; the fiducial parameters do not necessarily correspond with the parameters that have the largest share in the Top 
Set. Regularization terms are maximum entropy method (MEM), total variation (TV), total squared variation (TSV), the £, norm (£i), and the 


weighted-£, (£ n): and these are described in M 87* 2017 IV. 


S/N (Rogers et al. 1995; Blackburn et al. 2020). The distribution 
of images in the resulting Top Sets on the real M 87* data is a 
proxy for the uncertainties due to different imaging strategies 
and assumptions. For each method we also defined a “fiducial” 
image by converting the model pyx of each Top Set image to 
an effective blurring width, averaged this width across all syn- 
thetic data sets, and found the set of parameters that produced 
the minimum effective blurring width. A summary of the Top 
Set parameters for each method is shown in Table 2. 

Validation of the Top Set parameter combinations was then 
performed by imaging the remaining six geometric models as 
well as a GRMHD snapshot image and ensuring that the result- 
ing images closely match their ground truths (see Appendix G.1 
and Fig. G.3). We also found very few artifacts in the image 


reconstructions across the Top Set, as shown by stacking Top 
Set images for each synthetic data set as shown in Fig. G.4. The 
detailed imaging strategies of the three pipelines are summarized 
in Sects. 5.1.1—5.1.3. 

As mentioned above, we also utilized two methods, THEMIS 
(Broderick et al. 2020a) and Comrade (Tiede 2022), that fol- 
lowed a Bayesian posterior sampling approach to image recon- 
struction. To quantify the imaging uncertainty, THEMIS and 
Comrade use Bayesian inference and cast VLBI imaging as a 
Bayesian inverse problem. Since these Bayesian methods do not 
require training for selecting model hyperparameters, all syn- 
thetic data models are treated as validation data sets. In addition 
to producing a best-fit image to the relevant visibility data, these 
approaches produce an image family that captures the image 
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Table 3. Initial DIFMAP geometric models. 


cres180 ring dblsrc disk cres-980 cres8 cres90 edisk point+disk point+tedisk GRMHD HOPS CASA 
April21 Bandl R-48  R-44 R-40 R-44  R-44 R-44  R-44  D-56 P P R-52 R-48 R-48 
Band2  R-48  R-44 D-68 D-72  D-68 R-48  D-68  D-56 P P R-52 R-48  R-48 
Band3  R-84  R-48 D-40 D-68  D-84 D-84  D-84  R-40 P P P R-48  R-48 
Band4 R-48 R-48 R-40 D-68  R-48 R-48  R-48  R-40 P P P R-48  R-48 
April25 Bandl | R-40  R-40 P D-72  D-56 D-56 D-56  D-56 P P R-52 D-56 R-48 
Band2 D-64 D-60 P D-72  D-56 D-60  D-64  D-56 P P R-48 D-56 D-56 
Band3 R-44 R-48 P D-68 | R-44 R-44  R-44  D-56 P P R-36 R-40  D-56 
Band4  R-44  R-44 P D-68 | R-44 R-44  R-44  R-36 P P R-36 R-40  D-56 


Notes. These models are chose based on closure phase fitting for the 11 synthetic data and the real M 87* data calibrated by HOPS and CASA. 
We note that R, D, and P refer to a ring, disk, and point source model, respectively. The diameters of the ring and disk models in units of was are 


noted. 


uncertainty arising from the measurement uncertainty of the data 
given the assumptions of the model. This permits us to quantita- 
tively assess the significance of image structures, array calibra- 
tion quantities, and the relationship between these features. More 
detailed descriptions of these methods are given in Sects. 5.1.4 
and 5.1.5. 


5.1.1. DIFMAP 


The CLEAN algorithm (e.g., Hógbom 1974; Clark 1980) has 
been widely used for image reconstruction in radio interferomet- 
ric imaging. This algorithm uses an inverse modeling approach 
to derive a sparse reconstruction on the image domain using the 
interferometer's point-source response (i.e., the dirty beam). The 
CLEAN algorithm used in this paper assumes that the sky bright- 
ness distribution is a collection of point sources. The imaging 
process involves generating point sources (CLEAN components) 
at the location of the brightness peak in the dirty image (defined 
as the 2D Fourier transform of the measured visibilities) and 
iteratively subtracting them until a stopping criterion is reached. 
The final image is obtained by convolving the CLEAN com- 
ponents with a Gaussian restoring beam and adding the resid- 
ual image to represent the residual noise. Restrictions, known as 
CLEAN boxes or CLEAN windows, can be placed on the area 
of the map in which the CLEAN components are searched and 
are especially important for data with sparse (u, v) coverage. 

We conducted imaging with CLEAN using DIFMAP 
(Shepherd 1997, 2011) based on the imaging pipelines utilized 
for the 2017 EHT images of M87* (M87* 2017IV) and of 
Sgr A* (Sgr A* 2017 IIT). We perform the DIFMAP analysis using 
data averaged to 10s, and perform self-calibration using the 
same averaging time-scale. The presented pipeline has a few 
minor modifications compared to the 2017 M 87* pipeline. We 
omitted setting the estimated flux density expected from a base- 
line of zero length since there are actual data from intra-site 
baselines at very short lengths. Based on the gain error budget 
presented in Appendix E, we set acceptable limits on amplitude 
gain correction factors during the self-calibration process, which 
were within the range of 0.5—2.0, instead of the 0.83— 1.2 range 
used in the 2017 pipeline. We found large gain correction factors 
for the GLT, LMT, and PV at least for several scans. With the 
addition of the GLT as a new station in the array, the (u, v) cov- 
erage of the 2018 EHT data was improved, especially on inter- 
mediate baselines. This helped to suppress the side lobes in the 
synthesized beam, improving the CLEAN image reconstruction. 
The pipeline surveyed five parameters: the total assumed com- 
pact flux density, cleaning stopping condition, relative weight 


A79, page 12 of 63 


correction factor for ALMA in self-calibration, diameter of the 
CLEAN window, and the power-law scaling of the (u,v) den- 
sity weighting function. The same parameter ranges used for the 
2017 M 87* imaging were used, with a compact flux density of 
0.4 Jy to maintain consistency with RML imaging methods. As 
explained in M 87* 2017 IV, the amplitudes measured from LMT 
suffer from relatively poor a priori calibration and thus require 
self-calibration with the initial imaging result. 

We implemented phase-only self-calibration with geomet- 
ric models to mitigate the impact of a priori calibration uncer- 
tainties in EHT measurements during the imaging process. Two 
strategies were employed for selecting the initial model in the 
parameter survey using synthetic and real data sets. The first 
strategy involved employing a point source with a flux density 
of 1 Jy positioned at the origin of the dirty image, similar to the 
approach used in the 2017 M 87* imaging (M 87* 2017 IV). The 
second strategy aimed to choose a better geometric model for the 
initial phase-only self-calibration, following the methodology of 
the 2017 Sgr A* imaging (Sgr A* 2017 III). For the latter strategy, 
synthetic visibilities were generated for various geometric mod- 
els, including a Gaussian with 15 uas FWHM (representing an 
unresolved symmetric model), a uniform disk with sizes ranging 
from 56 to 84 uas in steps of 4 uas, and a uniform ring with sizes 
ranging from 36 to 68 uas (also in steps of 4 uas, without width). 
The best model was selected based on the closure phase nor- 
malized x2. Unlike the Sgr A* imaging, we did not incorporate 
further CLEAN reconstruction after determining the best initial 
model because it resulted in instability in the best initial model 
for the data presented in this paper. The initial models chosen 
for each data set can be found in Table 3. Our analysis revealed 
that selecting the geometric model based on closure phase fit- 
ting outperformed the point source model strategy for synthetic 
data reconstruction, yielding more Top Set images. This outcome 
was anticipated, considering the more complex source geome- 
tries in synthetic and real M 87* data compared to a simple point 
source. Even though the point source model strategy produced 
fewer Top Set images, the images themselves look similar to the 
Top Set images produced using the closure phase fitting strategy. 
We present the results from closure phase fitting in the main text 
and include the point source modeling strategy in Appendix H. 

The DIFMAP pipelines presented in this paper have been 
updated compared to those used in previous EHT imaging analy- 
ses (M 87* 2017 IV; Sgr A* 2017 III). We performed a backward 
compatibility test for these pipelines on the 2017 M 87* data 
(April 11, low-band). We found that the resulting images agree 
with those presented in M 87* 2017IV. We present the fiducial 
images of the M 87* data and the synthetic data in Appendix I. 
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5.1.2. eht- imaging 


The Python package eht- imaging (Chael et al. 2018, 2019a,b) 
is an RML-based VLBI imaging software capable of producing 
images by placing different relative weights on the fits to clo- 
sure quantities and complex visibilities. The results in this paper 
were produced by performing a parameter survey using the 2017 
M 87* imaging pipeline?. The surveyed parameters consist of 
the total assumed compact flux density, the fractional system- 
atic error on the measured visibilities, the FWHM of the circular 
Gaussian used as the initial and prior image, and weights for four 
of the regularizers. See Table 2 for a summary of the surveyed 
parameter space and Appendix A in M 87* 2017 IV for detailed 
definitions of each regularizer. The range of surveyed compact 
flux density values were chosen based on the pre-imaging con- 
straints outlined in Sect. 4, while the range of surveyed regular- 
izer weights were chosen based on experience from values sur- 
veyed in M 87* 2017 IV. All images were reconstructed with a 
128 uas FOV and a 64 x 64 pixel grid. 

The imaging pipeline starts by loading and coherently scan- 
averaging the data. Then the correlated flux densities at intra- 
site baselines were rescaled by the given compact flux density 
to remove the contributions from unresolved extended emis- 
sion outside the FOV. We added an additional fractional sys- 
tematic error term to the visibilities’ error budget to account for 
unknown amounts of non-closing errors in the data. As measure- 
ments taken by the LMT suffer from large gain uncertainties, we 
performed an initial amplitude-only self-calibration to the LMT 
data and adopted station based gain corrections for LMT. This 
self-calibration is performed with a circular Gaussian geometric 
model with FWHM of 60 pas and flux density of 0.6 Jy, chosen 
to fall in the center of the compact flux density limits derived 
in Sect. 4. Lastly, the visibility amplitudes were inverse tapered 
with a 5 uas FWHM circular Gaussian to enforce an angular res- 
olution limit on the final reconstructed image. 

After these pre-imaging calibration steps, the pipeline pro- 
ceeded with four iterations of imaging and self-calibration. The 
imaging was initialized with a circular Gaussian of FWHM and 
compact flux density specified by the given parameter combina- 
tion. The details of the self-calibration and the relative weights 
placed on fits to the various data products were modified between 
each iteration to reflect progressing amounts of confidence in the 
gain and phase solutions. The first two rounds of self-calibration 
were performed only on the phases while the last two rounds 
were performed on amplitudes and phases. For the relative data 
weights, we began the first round of imaging by placing unity 
weight on the closure quantities, a fifth of that on the visibility 
amplitudes, and no weight on complex visibilities. As we pro- 
gressed through iterations, we removed weight on the visibility 
amplitudes and allowed non-zero weight on complex visibilities. 
The ratio between weights placed on close quantities compared 
to complex visibilities decreased in later iterations as we con- 
verged on a phase solution. All self-calibration rounds were per- 
formed on the scan-averaged time intervals. 

Each iteration involves several attempts at producing an 
image to prevent the imaging function from getting stuck in a 
local minimum. Each attempt utilized the previous best-fit image 
blurred to the nominal array resolution as the initial image. At 
the end of all four iterations of imaging, we ensured consistency 
with the original data and limited the angular resolution of recon- 
structed features by convolving the final image with the same 


? https://github.com/eventhorizontelescope/ 
2019-D01-02/tree/master/eht-imaging 


5 was Gaussian used for inverse tapering in the pre-imaging cal- 
ibration step. 


5.1.3. SMILI 


The imaging pipeline SMILI (Akiyama etal. 2017a,b) is a 
Python and FORTRAN RML-based imaging software. Similar to 
the survey conducted with eht-imaging, the SMILI parame- 
ter survey allows variation of the following parameters: the total 
assumed compact flux density, the FWHM of the circular Gaus- 
sian used for the prior image, the fractional systematic error on 
the measured visibilities, and the weights of three regularizers. 
See Table 2 and Appendix A in M87* 2017IV for details of 
the SMILI parameter ranges and descriptions of the regularizers, 
respectively. All images are reconstructed with a 128 tas FOV 
and a 64 x 64 pixel grid. 

Before imaging, the script coherently scan-averaged the vis- 
ibilities. Pre-calibration of the LMT gains and the addition of 
non-closing systematic errors were performed as described in 
the second paragraph of Sect. 5.1.2. We performed four imaging 
cycles for each self-calibration stage, self-calibrating only to the 
final reconstructed images in each of the three cycles. Recon- 
structions at each stage were initialized with a circular Gaussian 
of FWHM = 20 yas containing the compact flux density of the 
given parameter combination. Fractional uncertainties of 50%, 
30%, and 5% were added in quadrature to the error of visibil- 
ity amplitudes on LMT, GLT, and baselines to other antennas, 
respectively, to account for residual errors in the amplitude cali- 
bration. After the first imaging attempt in each stage, subsequent 
initializations used the previously obtained image re-centered 
to the image center of mass and blurred with a 20 was FWHM 
circular Gaussian. Each imaging cycle performed 1000 itera- 
tions of the L-BFGS-B algorithm (Byrd et al. 1995; Zhu et al. 
1997) used for the gradient-descent optimization in SMILI’s 
image solver. After the three imaging cycles, complex visibil- 
ities were self-calibrated with their output image. In the first 
two self-calibration stages, the imaging step used closure data 
(closure amplitudes and phases) and visibility amplitudes. The 
imaging used closure data and complex visibilities for the final 
two self-calibration stages. Similar to eht-imaging, SMILI 
also performed all rounds of self-calibration on scan-averaged 
timescales. 


5.1.4. THEMIS 


We employ THEMIS, a generic modular parameter estima- 
tion framework specifically developed to compare parameterized 
models with VLBI data produced by the EHT (Broderick et al. 
2020a). The basic THEMIS image model aims to reproduce 
the on-sky brightness distribution and the time-dependent com- 
plex station gains. For reproducing generic brightness distri- 
butions, THEMIS utilizes an adaptive splined raster model 
defined by a set of brightness control points arranged on an 
adjustable rectilinear grid (Broderick et al. 2020b), which we 
call a “Themage” (THEMIS image, see Appendix J for details). 
We used a raster resolution of 5 x 5, which was chosen after 
a survey of other resolutions ranked by the Bayesian evidence 
and described in more detail in Appendix K. Since there was a 
significant difference in the flux density between the intra-site 
baselines and the shortest inter-site baselines, we also included a 
large scale asymmetric Gaussian component to approximate the 
flux density contribution from unresolved emission outside of 
the FOV. 
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The THEMIS image model’s raster and Gaussian com- 
ponents are constructed in the visibility domain and directly 
fit against the scan averaged full complex visibilities gener- 
ated from the standard 10s averaged network calibrated data 
described in Sect. 3.3, to which we added an additional 1% frac- 
tional systematic error in quadrature to the thermal uncertain- 
ties. When fitting complex visibilities, thermal errors are strictly 
Gaussian and improve the smoothness of the likelihood sur- 
face. When we reconstructed the scan-averaged station gains, 
we imposed Gaussian priors on the logarithmic gain amplitudes 
and flat priors on the gain phases. Informed by the analysis of 
crossing baseline tracks and the EHT station flux density calibra- 
tion parameters (Table A.1), for the April 21 data, we imposed 
prior widths of 10% for all stations except LMT and GLT, which 
used 30% and 100% prior widths respectively. We used the same 
gain priors for the April 25 data except for permitting a 100% 
gain prior for PV, since it was noted that PV exhibited large 
amplitude fluctuations after UTC 02:00 due to poor weather (see 
Appendix D). 


5.1.5. Comrade 


Comrade is a Bayesian differentiable modular modeling frame- 
work written in the Julia programming language for use with 
VLBI (Tiede 2022). In this work, we applied non-parametric 
modeling and include a rasterized image model similar to the 
one described in Broderick et al. (2020b). We first scan-averaged 
the standard 10s averaged network calibrated data described in 
Sect. 3.3. We added an additional 1% fractional systematic error 
in quadrature to the visibility thermal error and then extracted 
the closure phases and visibility amplitudes. We removed clo- 
sure phases with S/N < 3 and any closure phases with intra-site 
baselines. 

Our model for all image reconstructions consisted of a ras- 
terized image model (of dimensions N, x N,), a Gaussian of 
FWHM = 1 mas, and scan-averaged station gains. The Gaussian 
was used to model the amplitudes of short intra-site (ALMA-— 
APEX, JCMT-SMA) baselines. In the image model, the grid of 
raster points was convolved with a third-order basis-spline (B- 
spline) kernel (Eq. (L.2)) to generate flux densities that smoothly 
varied in all directions. Comrade's raster model is similar to 
THEMIS in Sect. 5.1.4, except that we used the B-spline ker- 
nel. A detailed description of the image model is given in 
Appendix L. 

The hyperparameters of the model are the FOV and raster 
size (N.XN,). For April 21 bands 3 and 4, we used a 12x12 pixel 
raster with a 7.5 uas pixel size, which was ~1/3 of the nominal 
resolution and FOV of 90 uas. This FOV was chosen by checking 
the visibility-amplitude residuals to incorporate all the flux den- 
sity within the FOV. If the FOV is too small («65 uas) then the 
residuals are large. If the FOV is too large (>90 uas), we do not 
have information on those scales from the visibilities and we get 
poor reconstructions for the synthetic data sets. For bands 1 and 
2, we shrunk the FOV to 75 was due to the lack of intermediate 
baselines from GLT; maintaining the same pixel size resulted in 
a 10x 10 raster. The hyperparameters are the same for M 87* and 
the synthetic data. The hyperparameters were changed depend- 
ing on (u, v) coverage for different bands and days, and the spe- 
cific details are given in Table 4. 

We formed the visibility amplitude likelihood and closure 
phase likelihood as described in Appendix F of Sgr A* 2017IV. 
We used a uniform distribution 44(0.0, 1.5) Jy for the prior on the 
total flux density f and U(0.0, 1.0) for the fraction of the total 
flux density f; for the flux density of the large-scale Gaussian. 


A79, page 14 of 63 


Table 4. Comrade raster model hyperparameters. 


Bands 1 and 2 Bands 3 and 4 


Day FOV (uas) N,xN, FOV (uas) N, x Ny 
April 21 75.0 10 x 10 90.0 12x 12 
April 25 67.5 9x9 67.5 9x9 


For the raster pixel fluxes, we used a symmetric Dirichlet distri- 
bution on the simplex (see Eq. (L.5)). For the station gain priors, 
we used a normal distribution for the log-gain amplitudes for 
each station. The prior width of station gain priors was similar 
to THEMIS in Sect. 5.1.4. A detailed description of the prior 
distributions is given in Appendix L. 

Finally, the un-normalized posterior was formed by taking a 
product of the prior and the likelihood distributions. To sample 
from the posterior, we used a two-stage strategy. First, to find a 
reasonable starting location, we adopted the L-BFGS optimizer 
(Liu & Nocedal 1989; Mogensen & Riseth 2018), running it 20 
times, and initializing it from random draws from the prior. We 
then selected the location that optimized the log joint probability 
density at our starting location for the Hamiltonian Monte Carlo 
No-U-Turn Sampler (Hoffman & Gelman 2014; Xu et al. 2020). 
We ran the sampler for 12000 steps; the first 10000 are adap- 
tation steps. To check for Markov chain Monte Carlo (MCMC) 
convergence, we measured the split-R (the average variance of 
draws in one chain compared to the variance across all chains) 
and the effective sample size. After sampling, we calculated 
statistics of the posterior, such as the mean and standard devi- 
ation of the images and station gains. Since the phase center for 
EHT-like data is unconstrained, we used the image centroid to 
align the images. 


5.2. Presentation of images 
5.2.1. Comparison between methods, bands, and days 


Figure 7 shows representative images of M 87* produced with 
each of the five imaging methods using data from both calibra- 
tion pipelines. Images are shown for all four bands on April 21 
and April 25. For the DIFMAP and RML methods, we present 
fiducial images. For THEMIS and Comrade, we display a ran- 
dom sample from the posterior. 

We first discuss images from the band 3 data on April 21, 
which represent (together with band 4) the best (u,v) cover- 
age and the most stable imaging results. On April 21, DIFMAP, 
eht-imaging, and SMILI could all produce a non-zero num- 
ber of Top Set images for all bands. A visual inspection shows 
that all images display similar characteristics, including diam- 
eter, a central flux depression, and a brightness asymmetry in 
roughly similar positions (see Appendix G.2). Apparent differ- 
ences in detailed structure between methods can be attributed to 
differences in the effective resolutions of the imaging pipelines. 
For example, a 20 was deconvolution beam was used for DIFMAP 
imaging, so DIFMAP images tend to have a larger ring width 
and weaker central depression. In addition to good agreement 
in structure among the image reconstructions, we also find good 
agreement in the reconstruction of the time-dependent station 
gains for both the synthetic and real data (see Appendix G. 1). 

Figures 8 and 9 show the data visibility amplitudes and 
closure phases with the corresponding model visibility ampli- 
tudes and closure phases computed from each April 21 band 3 
image. These figures demonstrate that the images produced by 
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Fig. 7. Representative images recovered from the HOPS and CASA data with all five imaging pipelines for two observing days (April 21 and 
25). Each panel shows the fiducial image of the corresponding top set images for the DIFMAP, eht-imaging, and SMILI pipelines, and a random 
sample from the respective posterior for the THEMIS and Comrade pipelines. We do not have Top Sets for band 1 and band 2 from DIFMAP, 
eht-imaging, and SMILI pipelines on April 25. The dashed horizontal line in each block separates the DIFMAP and RML methods above from 
the Bayesian methods below. The circles in the DIFMAP panels represent an effective Gaussian blurring kernel of 20 was. The solid lines in the 
THEMIS and Comrade panels represent the size of the blurring kernel used to achieve the same effective resolution as the DIFMAP method. 


all pipelines fit the data well. When comparing the model clo- 
sure quantities of a Top Set image against the real data closure 
quantities, we used the mean squared standardized residual (7), 
normalized by the number of data points to be compared. For 
example, we define the normalized closure phase Xæ as: 


1) 


where N,, is the number of closure phase data points, We is the 
model closure phase, vc is the real data closure phase, and cy, is 


the closure phase standard deviation. Similarly, we can construct 
the normalized log closure amplitude x2 ca as: 
og CA 


1 |In Âc - In Acl 
5 ; 
O'inAc 


2 
Mog CA 7 M n (2) 
nAc 
where Nn ac is the number of log closure amplitude data points, 


Ác is the model closure amplitude, Ac is the real data clo- 
sure amplitude, and o? Ac is the log closure amplitude stan- 


dard deviation. Readers can refer to Sect. 2.1 of M 87* 2017 IV 
for additional details about constructing normalized closure y? 
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Fig. 8. Closure phases plotted as a function of GMST on three selected triangles from the April 21 band 3 observations (black points). The error 
bars on the data points denote the +1o uncertainties. The colored and dashed lines indicate the model closure phase curves from the fiducial 


images or posterior samples produced by the five imaging pipelines. 


values, which should not be confused for the formal definition 
of a reduced y. In Table 5, we show the normalized X? val- 
ues for the fiducial images and normalized X? statistics across 
the Top Sets for DIFMAP, eht-imaging, and SMILI. The fidu- 
cial images from the RML methods are consistent with the data 
roughly within the thermal noise, and the normalized X? values 
have little scatter across the Top Sets. We added 10% systematic 
uncertainty to the real data in the evaluation process for DIFMAP 
because the image generation process did not use closure quanti- 
ties (unlike the RML methods), and down-weighted the ALMA 
visibilites during self-calibration, resulting in comparable nor- 
malized y? values to the RML methods. 

For the Bayesian imaging methods, the fit quality is assessed 
through comparisons of the complex visibility reduced y? for 
THEMIS, the visibility amplitude and closure phase reduced 
x’ for Comrade, log-likelihoods, the logarithm of the Bayesian 
evidence (log(Z)), and distribution of the residuals. For this fit 
quality assessment, the maximum likelihood model from the 
sampling chains were used for THEMIS and maximum a pos- 
teriori model for Comrade. We present the relevant reduced y? 
values in Table 6. 

We calculate the reduced x? by dividing the standard x? by 
the sum of the number of independent model and gain parame- 
ters subtracted from the number of data points, or for the case of 
the complex visibility reduced y? (Red. ANS 


Red. Ny = Ny — Nmodel DOF — Nain DOF; (3) 
1 |i - vp 

Red.y2 = ————- ) ——— 4 

CXV = Deed. Ny) 2, ol P 


where Ny is the number of complex visibility data points, 
Nmodei DOF i$ the number of model degrees of freedom, Ni por 
is the number of gain degrees of freedom (we multiply the final 
reduced number by 2 since we fit both the real and imaginary 
components), V is the model complex visibility, V is the data 
complex visibility, and oy is the complex visibility standard 
deviation. One can construct similar reduced y? quantities for 
the visibility amplitudes (Red. xà) and closure phases (Red. 
Xin) using the appropriate model and data products, and prop- 
erly counting the number of data points, model, and gain degrees 
of freedom to construct the reduced denominator. 

The THEMIS and Comrade raster reconstructions produce 
low reduced y? across all bands and both days. Since the 
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Comrade raster used a substantially larger raster grid than 
THEMIS, we expected the Comrade fits to produce reduced y? 
well below unity. In the THEMIS fits, we find that band 1 recon- 
structions feature systematically worse performance. This may 
be indicative of additional systematic errors in those data sets. 

Details of the observations contribute to the differences 
between images from different bands and days. The improved 
(u, v) coverage in bands 3 and 4, given by the participation of 
the GLT, allows for improved reconstructions of M 87* images. 
The GLT is especially important in probing the null point near 
3.4 GA. This is proven by the increased number of Top Set 
images for the trained methods and the cleaner reconstruction 
of the ring-like morphology compared to bands 1 and 2. 

On the April 25 data, DIFMAP, eht-imaging, and SMILI 
struggled to produce a significant number of Top Set images — 
none of the methods could produce Top Set images for bands 1 
and 2. The Bayesian methods also struggle to produce a ring-like 
morphology for data taken on this day. This performance issue 
is mainly due to a lack of data from LMT on this day. The LMT- 
SMT baseline provided the only probe of the visibility structure 
around 1 G4; the lack of this baseline hampered imaging, espe- 
cially the capacity to select a preferred compact flux density. 


5.2.2. Compact flux density for RML and CLEAN 


The M 87* fiducial images were reconstructed with a total com- 
pact flux of 0.4Jy for DIFMAP, 0.5Jy for eht-imaging, and 
0.6Jy for SMILI; all of these values fell within the range of 
allowed compact flux densities derived from the pre-imaging 
constraints (Sect. 4). The RML methods only modeled the com- 
pact emission (since they rescale the intra-site baselines), so it is 
necessary to confirm that the image structure would not signifi- 
cantly change within the range of compact flux densities allowed 
by the pre-imaging constraints. Though each pipeline surveyed 
over several different assumptions on the total compact flux den- 
sity (see Table 2 for a summary), we found no strong prefer- 
ence for a particular value. As a complementary check we gen- 
erated images without intra-site visibilities and verified that the 
image structure is not significantly impacted. See Sect. 7.3 for 
a related discussion of the compact flux density estimates from 
the Bayesian methods, and Appendix M for supplementary syn- 
thetic data validation tests at different compact flux densities for 
the RML and CLEAN methods. 
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Fig. 9. Visibility amplitudes of band 3 data on April 21 as a function of baseline length compared with corresponding gain-calibrated visibility 
amplitudes from the five representative image models from each pipeline. The fiducial images are used for DIFMAP, eht-imaging, and SMILI, and 
the maximum likelihood model from the sampling chains are used for THEMIS and maximum a posteriori model for Comrade. The gray points 
represent the data used by each method after flagging, but before individual self-calibration. eht-imaging scales the intra-site baselines to 0.5 Jy, 
and SMILI scales the intra-site baselines to 0.6 Jy, and both pre-calibrate the LMT-SMT baselines before fitting. DIFMAP uses 10 s averaged data, 
and all other methods use scan-averaged data. THEMIS and Comrade apply 1% fractional systematic noise to all baselines, and eht-imaging 
and SMILI apply 1% and 0% fractional systematic noise respectively to all baselines, seen as minor differences in the gray points. The colored 
points represent the image model visibilities from each method, with station gains derived from each method's internal self-calibration procedure 
applied to the image model visibilities. Below each visibility amplitude figure are the normalized residuals for each image, which is the difference 
between the gray and the colored points, divided by the uncertainty of the gray data points. 
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Table 5. Fiducual image and Top Set closure quantity normalized y? 
values. 


DIFMAP eht-imaging SMILI 
Fiducia: — yt, 0.94 1.09 123 
Xing ca 046 1.01 1.23 
Top Set: y2p 0.91+0.33 1.24+0.29 1292027 
Xi, ca 0.89+0.39 1.16+0.28 1.37 + 0.30 


Notes. Closure quantity normalized x? values for the fiducial M 87* 
images and normalized y? statistics (mean and standard deviation) for 
Top Set images from the April 21 band 3 observations. We show the 
results assuming systematic uncertainties of 10% for DIFMAP and 0% 
for the RML methods. 


The difference of the flux densities between F;y; and Fepct 
implies there is some emission outside the compact region. 
However, its structure (e.g., scale and direction) is uncon- 
strained by the data (Sect. 4). The sub-milliarcsecond scale 
jet is not detected by the 2017 and 2018 EHT data, but is 
clearly seen in the 3mm images (Hada et al. 2016; Kim et al. 
2018a; Lu et al. 2023). This sub-milliarcsecond scale jet is pre- 
sumably resolved out due to lack of visibility and closure 
information from short enough (u,v) spacing (e.g., x1G4A). 
The contribution from the sub-milliarcsecond scale jet becomes 
less important at (u,v) distance 21 GA in 3mm data (Lu et al. 
2023). Assuming the jet emission is optically thin synchrotron 
radiation (EHT MWL Science Working Group 2021), the sub- 
milliarcsecond jet can be less bright at higher frequencies, so we 
expect that it would be challenging for the EHT to detect without 
much better (u, v) coverage at short baselines. 

As we did in M 87* 2017IV, the geometric models in the 
synthetic data sets included emission from a simple extended 
jet modeled out to ~2 mas. None of the image reconstructions 
on the synthetic data recover this extended component. Simi- 
larly, no method was able to recover diffuse or extended emis- 
sion displaced from the ring seen in the GMRHD synthetic data. 
While several studies suggest the presence of a small-scale dif- 
fuse structure in the 2017 EHT data reconstructions which aligns 
with the limb-brightened jet observed at 3 mm (Arras et al. 2022; 
Carilli & Thyagarajan 2022; Broderick et al. 2022a), no con- 
straints on larger-scale (70.2 mas) jet emission are supported by 
the data. 


5.2.3. Image statistics 


The distributions of both the Top Set and posterior images 
help us understand the image uncertainties associated with each 
method. Figures 10-12 show the image- and visibility-domain 
uncertainties associated with the image sets for eht-imaging, 
THEMIS, and Comrade, respectively. 

The uncertainties shown for eht-imaging are a proxy for 
the uncertainties in choosing the regularizer weights and hyper- 
parameters. Similar to the corresponding figure in M 87* 2017 IV, 
we find that the high image uncertainties correspond to locations 
of high brightness temperature and visibility-domain uncertain- 
ties primarily due to gaps in (u, v) coverage. We see in the image 
domain figure of normalized standard deviation that small con- 
centrations of uncertainties exist at various locations along the 
ring. However, they are less pronounced than the "knots" in 2017. 
The very small amount of uncertainty in the central flux depres- 
sion indicates the robustness of the ring-like feature, as it occurs 
in nearly every Top Set image. The normalized standard devi- 
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ation visibility domain image shows a similar concentration of 
uncertainty around 2 GA. The fractional standard deviation image 
shows a sharp increase in uncertainty at the boundary between 
(u, v) distances probed by the EHT versus those outside the max- 
imum probed distance. This indicates that the RML methods are 
not assigning high confidence to image features smaller than the 
minimum scale probed by the observations. 

Comparatively, the uncertainties derived from THEMIS and 
Comrade stem from statistical uncertainty of the image pos- 
terior space given the model assumptions. Both THEMIS and 
Comrade find that the ring-like feature is robustly recovered, 
finding it has a low fractional standard deviation. Additionally, 
both THEMIS and Comrade find no evidence for non-ring emis- 
sion. For THEMIS, this can be seen from the fact that the raster 
is concentrated around the ring. Substantial non-ring emission 
would cause the THEMIS raster to expand to cover the region, 
as seen in Fig. G.3. Comrade's larger assumed FOV allows for a 
direct quantitative estimate of the significance of non-ring emis- 
sion and finds no significant (>3c) detection of non-ring emis- 
sion. Moving to the visibility domain, we see that the uncer- 
tainties reported by THEMIS and Comrade are qualitatively 
different and generally find that THEMIS’ uncertainty is much 
smaller than Comrade's. This disparity reflects the differences 
in each method's model specification. For short (u, v) distances, 
the lower standard deviation for THEMIS is due to its sparsity- 
inducing prior and fitted raster dimension, which tends to pro- 
duce smaller image FOVs than Comrade, reducing the uncertainty 
on scales <3 GA. Additionally, since THEMIS uses a relatively 
large pixel size (212 uas), itis unable to produce significant ampli- 
tudes on scales 210 GA, which drastically decreases its measured 
visibility-domain uncertainty compared to Comrade. 


5.8. Image domain feature extraction (IDFE) 


While the imaging methods described above produce excellent 
fits to the visibility data and produce a consistent ring-like struc- 
ture, the specific choice of regularizers and weights, masks, and 
parameter combinations used to produce these image reconstruc- 
tions are generally agnostic to the image morphology and the 
hyperparameters may not directly relate to physical quantities. 
To produce measurements of ring diameter, width, flux asym- 
metry, depth of the central depression, etc., we pass the recon- 
structed images through feature extraction algorithms similar 
to the analyses featured in M 87* 2017 IV and Sgr A* 2017 III. 
The entire Top Set is used for the CLEAN and RML meth- 
ods to obtain these measurements, while 500 random samples 
drawn from the respective posteriors are used for THEMIS 
and Comrade. Negative-valued pixels in THEMIS images are 
replaced with zeros in line with the previous Sgr A* analysis 
(Sgr A* 2017 HI; Sgr A* 2017 IV). 

We use two image domain feature extraction tools, REx 
and VIDA. The Ring Extractor (REx) is available as part of 
eht-imaging and is described in detail in M 87* 2017 IV. The 
second IDFE tool, Variational Image-Domain Analysis (VIDA, 
Tiede et al. 2022a), uses parameterized templates to approximate 
an image and adjusts the parameter values until a specified cost 
function, in the form of a probability divergence that provides a 
distance metric between the image and template, is minimized. 
We use an updated version of the SymCosineRingwFloor 
template used in Sgr A* 2017IV to match the mF-ring model 
described later in Sect. 6.2. We set m = 2 (the azimuthal 
brightness mode) for consistency with the geometric modeling 
analyses (see Sect. 6.1). While REx is intended to work only on 
ring-like images, VIDA can support a variety of templates, which 
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Table 6. Reduced x? quantities for the THEMIS and Comrade raster models. 


THEMIS raster Comrade raster 
Red. P Red. Noe Red. Xe Red. P m 
April21  Bandl 1.41 0.64 0.44 0.50 
Band 2 1.1 0.80 0.35 0.47 
Band 3 1.14 0.94 0.56 0.69 
Band 4 1.19 0.92 0.49 0.63 
April 25 Band 1 2.20 0.93 0.32 0.43 
Band 2 1.71 0.80 0.28 0.35 
Band 3 1.17 0.82 0.37 0.37 
Band 4 0.96 0.84 0.52 0.36 


Notes. For THEMIS, we report the complex visibilty reduced y? for the best fit model to the HOPS data. For Comrade, we report the closure phase 
reduced y^, the visibility amplitude reduced 7, and the combined closure phase and visibility amplitude reduced y? for the best fit model to the 
HOPS data. 
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Fig. 10. Visualization of image statistics calculated using the Top Set images from the eht-imaging pipeline for observations taken on April 21 
band 3. We emphasize that these images do not represent the posterior probability space for the reconstructions. Each image reconstructed using 
eht-imaging is the maximum a posteriori (MAP) image for a given parameter set. Thus, the statistics shown represent uncertainties that arise 
from different choices of regularizer weights, not from an exploration of posterior space. The top row shows top statistics in the image domain 
while the bottom row shows the visibility domain. Overlaid on the visibility domain panels is the (u, v) coverage for the April 21 observation. 
From left to right, we present the mean image; the standard deviation; the normalized standard deviation, calculated by rescaling each image to the 
flux of the mean image; and the fractional standard deviation, calculated by dividing the standard deviation by the mean. The fractional standard 
deviation panel has been clipped to a maximum value of 1. Portions of the image exhibit large fractional standard deviations due to pixel values 
very close to zero in the mean image. In the top row, image contours are drawn at 10%, 20%, 40%, and 80% of the peak values from the mean 
image. In the bottom row, the gray contours represent 0.1%, 1%, and 10% of the peak while the black contours represent 10 and 100 mJy (left 
three panels) and 0.1 (right most panel). The complex visibilities are calculated by taking a Fourier transform of the images and then calculating 
the mean and standard deviation. The absolute value of the mean and standard deviation of the complex visibilities are used to calculate visibility 
amplitudes. 


can be used to evaluate the success of our imaging methods by 
quantitatively reproducing the non-ring synthetic data sets such 
as the disk and the double source. 

Figure 13 shows the ring parameters (diameter d, width w, 
and position angle) measured with REx and VIDA for all bands, 
days, and imaging pipelines for HOPS and CASA data recon- 
structions. On April 21, all bands and pipelines show ring- 


like structures with a roughly consistent diameter of ~42 uas 
as expected from the image morphology discussed in Sect. 5.2. 
Position angle values are around 215? and are consistent among 
all bands and pipelines for April 21. Width measurements have 
relatively large differences among the imaging pipelines, espe- 
cially between DIFMAP and the other methods. This is primarily 
due to the beam convolution effect in CLEAN imaging. Since 
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Fig. 11. Visualization of image uncertainties using images from the posterior of THEMIS pipeline for observations taken with band 3 on April 21. 


The contour lines shown are drawn as described in Fig. 10. 
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Fig. 12. Visualization of image uncertainties using images from the posterior of Comrade pipeline for observations taken with band 3 on April 21. 


The contour lines shown are the same as described in Fig. 10. 


measurements of the ring width are dependent on an individ- 
ual method’s ability to leverage super-resolution, it is challeng- 
ing to combine the width measurements across methods. While 
some methods trend toward slightly lower widths compared to 
their 2017 values (DIFMAP, eht-imaging, THEMIS), others are 
more consistent with the 2017 values (SMILI, Comrade). While 
the results from two independent IDFE methods, REx and VIDA, 
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generally agree, there do exist some discrepancies. These dis- 
crepancies between the IDFE methods primarily arise due to the 
difference in the functional responses to an image with double 
rings or otherwise corrupted ring-like structures. 

Tables 7—10 list the diameter, width, position angle, asym- 
metry, and fractional central brightness measured with REx and 
VIDA from all days, bands, a priori calibration, and imaging 
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Fig. 13. Ring characteristics from all bands, observational days, and imaging pipelines coming form HOPS and CASA data (colored points). 
Median values and 68% confidence intervals are shown for diameter and width. Circular mean and standard deviation values are shown for the 
position angle. Circle and triangle markers correspond to REx and VIDA respectively. Red (eht-imaging) and blue (DIFMAP) dashed lines and 
shaded regions show the ring parameters and 68% confidence intervals from the 2017 April 11 measurements. The vertical gray shaded region 
corresponds to the 2018 April 21 measurements, and the vertical un-shaded region corresponds to the 2018 April 25 measurements. 


pipelines. The brightness asymmetry A is ~0.2—0.4, indicat- 
ing an asymmetric ring-like structure. Moreover, the fractional 
central brightness fc is —0—0.5, indicating a clear central 
depression. 

Comparing the extracted features between April 21 and April 
25 is challenging due to the lack of Top Set images for DIFMAP, 
eht-imaging, and SMILI on April 25 bands 1 and 2 as well 
as the sparser (u,v) coverage on April 25 leading to higher 
uncertainty. For the available measurements, the features look 
consistent within their errors between days. Major discrepan- 
cies in the measurements, especially in the ring width, arise from 
the broken ring structures, as shown in Fig. 7. Ring parameters 
between HOPS and CASA agree well with each other, and we 
obtain a diameter of ~42 uas and position angle of ~215°. All 
parameter measurements except for the position angle are con- 
sistent between 2017 and 2018. The position angle has changed 
from ~180° in 2017 to ~215° in 2018. The interpretation of this 
change is discussed in Sect. 7. 

In summary, the 2018 EHT data is well described by a ring- 
like crescent structure. Where the data quality and coverage are 


good, we find a consistent diameter compared to the 2017 images 
and a slightly shifted position angle. Where the (u, v) coverage is 
poor, such as in bands 1 and 2 on April 25, we expect our recon- 
struction methods to struggle, and indeed we have much larger 
uncertainties in the image structure for those data sets. Now that 
we have determined the presence of a ring-like structure in the 
2018 data, in the next section we can compare and verify the 
measured image domain parameters by directly modeling those 
features with an assumed ring in the reconstructions. 


6. Direct modeling of physical parameters 


In addition to the image-domain feature extraction methods, 
we also attempted to model the physical ring parameters from 
the visibility domain information directly. We use a proce- 
dure similar to one described in M 87* 2017 VI where we con- 
struct a family of geometric models which closely match the 
observed M87* visibilities. This procedure can be seen as 
another form of image reconstruction with stronger assumptions 
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Table 7. Ring parameters for 2018 April 21 HOPS data. 


d (pas) w (uas) Position angle (°) A fe 
DIFMAP 
April 21 Bandi REx 3865331 27.2422 218.3 + 7.2 0.271%% 0.484015 
VIDA 40.8113  242:13 214.7 + 2.5 0.4340% — 0,52*007 
Band 2 383:23. 27,353 215.0 + 9.6 0.23008 — 9,493016 
40.813 242707 212.7 + 2.8 0.417955 — 0.51010 
Band 3 392043. 29192 217.7 + 13.8 0.23:000 — 0.46012 
39.715 2451$ 2152 + 2.9 0.307000 0.471913 
Band 4 39.1*1$ 25.6173 216.2 + 4.4 0.2900? — 0.407007 
39471] 23.713 213.6 + 22 0.35'000 0.427708 
eht-imaging 
April 21 Band! REx  412*1$ 140*H 218.4 + 8.9 0.203009 — 0,027002 
VIDA  415'0] 13.4#18 217.3 € 2.8 (28:00 9.01700 
Band 2 41.8189 13.545 214.6 + 10.0 0.277354 0.017581 
42.553 — 13.4722 213.2 + 3.0 0.297}; 0.017007 
Band 3 40.3:09 12.8712 220.7 + 24 0.2700? — 0.007005 
40.707 — 12.715 216.9 + 3.0 02700? — 0017007 
Band 4 40.705 12.23 2134 + 42 0.29733 — 0.007005 
413'80 12543 209.6 x 2.3 028'00 0.017007 
SMILI 
April 21 Band! REx  41.0*17 15.8535 220.8 + 9.7 0.24*000 — 9,03*01] 
VIDA  413*16 147333 217.6 + 3.0 0.21:009 — 0,037010 
Band 2 41329 15453 216.5 € 11.5 027009 00100 
41912. 14555 216.5 € 2.9 0.257353 — 0.027005 
Band 3 40.1413 14.723 220.4 + 7.0 0.277} 53 — 0.007000 
40.907 143732 2144 + 3.1 0.24'00 — 0,027008 
Band 4 39.8418 — 15.125 219.1 + 7.2 0.307}, — 0.007005 
40.8'10 — 14.8722 214.0 € 2.7 0.28001 — 0.037006 
THEMIS 
April 21 Band! REx  422*2 13.4715 220.8 + 6.0 0.27:003 — 0,027006 
VIDA  429*3 139719 2224 + 3.0 0.31*007 — 0,02*005 
Band 2 43.3415 13.243 220.9 x 3.5 0.28'009 — 0.007002 
43.3417 13.62 223.8 + 3.1 0.307000 — 0.017007 
Band 3 43.703. 11.0705 216.3 + 1.2 0.357}; 0.00754 
43.5'09 — 10.5705 216.9 + 2.5 0.36002 — 0.007000 
Band 4 434*0$ 11.844 212.6x 1.1 0.36001 — 0.007005 
43.659 uani 215.2414 0.40'09 — 0.01700) 
Comrade 
April 21 Bandl REx 41314  14.8*19 222.1 + 3.0 0.30*009 — 0.147006 
VIDA 42.1415 15.012 220.7 + 3.0 034'00$ 0.11005 
Band 2 42.9t18  145' 217.6 + 5.7 0.307000 — 0.117003 
42.6521 15.513 217.2 + 3.1 0.33000 — 0.097007 
Band 3 42.22 13.070? 2132 + 4.9 0.217953 — 0.08005 
42.5110 — 5.5102 2122 + 2.0 0.26000. — 0.04700) 
Band 4 43711. 14.2408 2124 + 6.0 0.23:009 — 0.08005 
44.3508 14405 213.6 + 1.4 0.29009 — 0.087005 


Notes. Median values and 6846 confidence intervals are shown for diameter, width, asymmetry, and fractional central brightness depression 
(d, w, A, fc). Circular mean and standard deviation values are shown for the position angle. 
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Table 8. Ring parameters for the 2018 April 25 HOPS data. 


d (uas) w (uas) Position angle (^) A fe 
DIFMAP 
April 25 Band 1 REx - - = = — 
VIDA - - - = = 
Band 2 - - - = = 
Band 3 36.808 22255 233.2 + 11.0 0.26007 — 0.237004 
37.9159 22653 231.0 € 2.0 0.25*} 97 — 0.27018 
Band 4 32.02 22,892 2143 + 41 0.33'00 — 0,347027 
34607;  222* 212.0 + 0.1 0.36000 — 0.4405 
eht-imaging 
April 25 Band 1 REx - - = = = 
VIDA - - - 2 7 
Band 2 - - - = - 
Band 3 369*4 — 16404 192.4 + 14.6 0.20'00? — 0.067002 
37T 45 — 19155 194.5 + 2.1 0.2800 0.19095 
Band 4 37.8108. 153109 207.3 + 14.9 0.28007 — 0.027001 
39.03. — 17.054 213.0 + 2.6 0.20*055 — 0.09*073 
SMILI 
April 25 Band 1 REx - - = = = 
VIDA - - = = = 
Band 2 - - - = = 
Band 3 38.5253 18.659 359.6 + 100.2. 0.32701 — 0.107026 
33.1774 — 3L350. 2715 + 2.3 0.25%} 56 0.679.453 
Band 4 38.25 — 16.8725 231.9 + 84.4 0.324012 0.037013 
38.4413, 22.173! 248.2 + 2.7 0.314919 — 0.39030 
THEMIS 
April 25 Bandi REx 38.0%} 15.15 271.7 + 59.4 0.307012  0.04+0.03 
VIDA 37.8432 16.8529 268.5 + 2.5 0.39%} — 0.187015 
Band 2 38755 — 17323 198.4 + 69.3 0.444575 — 0.04700 
34.737 — 18.6525 185.8 + 1.3 0.46750 0.197979 
Band 3 43.595 — 12.6500 218.2 + 4.5 0.31000 0.01460 
41655 B 217.4 + 2.5 0.33'00$ — 0.0100 
Band 4 398'1] 16253 219.3 + 15.3 0.2900% — 0.017001 
407326 — 173754 226.7 x 2.5 0.3600 — 0.08701! 
Comrade 
April 25 Bandl REx  432?]  152*53 2443 x 21.3 0.21000 — 0.24008 
VIDA 447553 — VL8'S 240.9 + 3.0 0.20*39 — 0.25700 
Band 2 425'31 149 224.0 + 28.6 0.20%} 93 0.22700 
4851. 17517 2222 + 3.0 0275 02205 
Band 3 42.6713 — 14.803 227.0 + 8.6 0.217305 0.1579 0 
434?! 161*1? 219.0 x 3.1 0:27:00. 0121007 
Band 4 421*13.— 14753 230.9 + 10.2 0.24005 — 0.167006 
41.5?) — I61*fi 224.8 + 3.0 0.2600; — 0.087005 


Notes. Median values and 68% confidence intervals are shown for diameter, width, asymmetry, and fractional central brightness depression 
(d, w, A, fc). Circular mean and standard deviation values are shown for the position angle. 
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Table 9. Ring parameters for 2018 April 21 CASA data. 


d uas) w (uas) Position angle (°) A fe 
DIFMAP 
April 21 Bandl REx 390532  29.4+17 220.1 + 20.5 0.257008 — (553012 
VIDA 40.3433 24.623 215.8 + 2.7 0.45*005 — 0,5900 
Band 2 34.457  29.8*09 208.2 + 41.4 0.16031 — 0.6570-18 
38.2718 26.071 205.7 + 2.1 0.334519 — 0.637057 
Band 3 38623. 2TLTU3 2184 + 23.3 0.22:000 — 0.49*010 
39.3432 25-739 216.0 + 2.8 0.32'0. 0.537017 
Band 4 37.9533 274028 216.1 + 10.3 0.22000 — 0.48013 
387335  265*H 213.2 + 2.8 0.29+001 — 0.50701] 
eht-imaging 
April 21 Bandl REx  405* 15.729 212.1 + 20.7 0.205005 — 0,0400 
VIDA 413309 15.6729 210.1 + 2.8 021700 — 0,051005 
Band 2 41.225) 144115 208.7 + 17.6 0.24'009 — 0.027005 
423*10 — 14,79 209.7 + 2.8 0.277}; 0.02700 
Band 3 39.703 13.3418 2154 x 23 0.24'001 — 0.027005 
404703 13.9418 212.8 + 2.6 0.24'00? — 0.02700 
Band 4 39.8759 134015 212.1 + 19.6 0.244} 7 0.017957 
40.9738 — 14.0723 212.9 + 24 0.25000 — 0.02700 
SMILI 
April 21 Bandl REx 40625 187553 2144 + 19.5 0.18199% 0.11297 
VIDA 422* 17.9723 218.6 + 3.0 0.24*005 — 0.163026 
Band 2 41.0479 17074 214.6 + 9.0 0.24'008 — 0.0870 
4174 16194 214.0 + 2.6 0.235009 — 0.067002 
Band 3 39.813 16.1725 219.7 € 7.7 0.227001 — 0.037008 
40.6505 — 16.055 213.8 + 3.1 0.21'000 — 0.057009 
Band 4 39.0413 15.621 218.5 + 13.5 0.227009 — 0.027008 
402*1 — 15.527 217.1 + 2.9 0.221901. — 0.057007 
THEMIS 
April 21 Bandl REx 380752 19.8704 198.3 + 18.0 0.19006 — 0.08*002 
VIDA 40.9733 22.531 211.1 € 2.5 0.235005 — 031045 
Band 2 429:20 14.102 216.9 + 5.6 0.25:003 — 0.0009] 
439*14 — 14.270? 2232 + 2.7 0.294021 9,040.05 
Band 3 43.5504 = 13,0704 213.1 1.5 0.321991 — 0.007000 
44.593 — 11,6503 219.0 + 1.5 0.505099 — 0.06700? 
Band 4 42.8:03 13.730; 213.5 + 1.3 0.314001 9,0070-00 
43513 13142 220.8 + 1.4 0.34*016 — 9,0270-06 
Comrade 
April 21 Band] REx  42.8:15 14.672 222.1 + 4.7 0.20004 — 0.167006 
VIDA 42.8413 17273 217.7 € 2.9 0.325004 — 0.13700 
Band 2 430*17 — 148*15 222.8 + 13.5 0.24*005 — 0.15*006 
45.0721 14,612 2214 + 3.0 0.36009 — 0.19*005 
Band 3 418113 12.9708 210.1 7.6 0.19009 — 0.08700 
42.510 — 16.1702 210.7 1.1 0.235009 — 0,05*002 
Band 4 424*16 — 14.0702 213.2 + 9.0 0.161993 — 0.107002 
43.513 — 16.805 2114 € 2.5 0.21008 — 0.06700 


Notes. Median values and 68% confidence intervals are shown for diameter, width, asymmetry, and fractional central brightness depression 
(d, w, A, fc). Circular mean and standard deviation values are shown for the position angle. 
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Table 10. Ring parameters for 2018 April 25 CASA data. 


d (uas) w (uas) Position angle (°) A fe 
DIFMAP 
April 25 Band 1 REx - - = = = 
VIDA = - = = = 
Band 2 - - = = = 
Band 3 35.9799 23.153 239.9 x 15.7 0.169005 — 0.1905 
38.44) — 23.3553 240.5 + 2.3 0.18%} 39 — 0.39013 
Band 4 35112 23.812 246.8 + 9.6 0.227003 0.3479 08 
36.07}4 — 24.6723 244.7 + 3.2 0.23°008 — 0.447008 
eht-imaging 
April 25 Band 1 REx - - = = = 
VIDA = = = = = 
Band 2 - € - = = 
Band 3 36559. q7719 2120 € 71 0.18700 — 0.05700 
3715'$ 20.3435 206.1 + 3.1 0.194303 — 0.22006 
Band 4 315349 3679 212.0 + 15.4 0.277003 — 0.047002 
38.2405 — 19.053 213.3 + 3.0 0.26¢ 351 oarn 
SMILI 
April 25 Band 1 REx - - = = = 
VIDA = - = = - 
Band 2 - e: = = - 
Band 3 39.5713 182739 261.4 + 242 0.30:006 — 0.06017 
39313 2282 257.8 + 2.2 0.37500 0.3419% 
Band 4 39.55 17.6438 244.4 + 27.4 0.257007 — 0.04701 
39.7713 — 20.819 234.8 + 3.0 0.3253] — 028015 
THEMIS 
April 25 Bandl REx  376'À! 162790 218.7 x 88.9 0.375021 — 0.067010 
VIDA 35.825 17.064 2202 + 2.1 0.42*008 023707 
Band 2 39.035 — 18.5770 197.2 + 53.5 0.43+020 0.047010 
33.5783 20.51 190.6 + 1.6 0.437993 — 03I'015 
Band 3 433*M. 13295 215.5 + 14.5 0.2670 — 0.00002 
43.5515. 13.1729 214.8 + 2.8 0.30709 — 0.0100 
Band 4 39.926 — 18.6523 244.6 + 65.1 0.19005 — 0.0805 
41.0:81 — 22.9555 242.5 + 2.9 0.25'009 — 03001) 
Comrade 
April 25 Bandl REx  425*4 — 153*H 250.5 + 36.0 0.7:007 0,261010 
VIDA 45055 17.6513 242.9 x 3.0 0.2790 — 028005 
Band 2 42156 i5] 219.7 + 34.8 0.197006 — 0,25100 
43.5533 — 18072) 220.3 + 2.9 0.2030 — 027010 
Band 3 428/53 1467] 219.1 + 30.4 0.16001. 0.2270 
44755  162*3 216.8 + 3.0 0.25*0]8 — (351015 
Band 4 45075 142] 218.5 + 12.0 0.18005 — 0.197009 
46.6735 — 168*11 215.5 + 3.0 0.39%} 13 — 0267005 


Notes. Median values and 6846 confidence intervals are shown for diameter, width, asymmetry, and fractional central brightness depression 
(d, w, A, fc). Circular mean and standard deviation values are shown for the position angle. 
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Fig. 14. Comparison of A log(Z) for a series of geometric models ordered by model complexity. The Bayesian evidence for each model is evaluated 
with data generated from closure amplitude and phases. The number of parameters needed to define each model is given in parentheses. Circle 
markers denote ring-like models and hourglass markers denote other models. colors denote models with similar construction. The right panel 


shows a zoom in of the gray shaded region of the left panel. 


on the basic image structure. First, in Sect. 6.1, we assess 
and quantify the extent to which ring-like features are pre- 
ferred by the data by comparing the Bayesian evidence of 
both ring-like and non-ring visibility-domain geometric mod- 
els. We then describe our fiducial models for direct physical 
parameter extraction and analyze their inferred ring features 
in Sect. 6.2. Finally, in Sect. 6.3, we check for consistency in 
the ring parameters across models and compare the results to 
M87* 2017 VI. 


6.1. Geometric modeling evidence exploration 


We estimated the preference for ring-like models by compar- 
ing geometric models of varying degrees of complexity to the 
2018 M87* EHT data. This is very similar to the procedures 
found in M 87* 2017 VI and Sgr A* 2017IV to evaluate the per- 
formance of different models against the EHT data. In this 
instance, we used different numbers and combinations of Gaus- 
sians, disks, rings, crescents, and m-rings (see Johnson et al. 
2020 and Sect. 4.3 in Sgr A* 2017 IV) as our model library. We 
evaluate each model based on the difference in the log Bayesian 
evidence, A log(Z). 

Our Bayesian evidence analysis was performed on data gen- 
erated from closure quantities constructed with band 3 observa- 
tions made on April 21 after applying an S/N > 3 cut. This 
means that these reconstructions are not susceptible to station- 
based corruption effects. The Bayesian evidences and posteriors 
for each model were evaluated using Comrade with a Nested 
Sampling (NS; Skilling 2006) posterior reconstruction scheme 
through dynesty (Speagle 2020). 

The resulting Bayesian evidence of each model is shown 
in Fig. 14, where it is evident that ring-like models provide a 
much better fit to the data when compared to non-ring models. 


A79, page 26 of 63 


We select the best-performing crescent and m-ring models; the 
THEMIS xs-ringauss model from M 87* 2017 VI (11 parame- 
ters) and a stretched m-ring of order 2 (8 parameters). We use 
these two models as the basis for defining fiducial models for 
direct feature extraction. 


6.2. Direct geometric modeling 
6.2.1. THEMIS direct modeling methods 


We expect that if the object we see at the center of M87 is a 
black hole (or sufficiently similar object), and the material sur- 
rounding it is optically thin, it should have a ring-like image fea- 
ture associated with half-integer orbits (n) of photons around the 
central gravitational object. Thus, the n = 0 “ring” corresponds 
to direct emission from around the black hole, the n = 1 ring cor- 
responds to photons that complete a half-orbit around the back 
of the black hole to impact our line of sight, and so on. This fea- 
ture is consistently seen in analytic, semi-analytic, and GRMHD 
simulations of optically thin accretion flows around black holes, 
and directly modeling this feature can, in principle, provide 
strong constraints on space-time properties (Johnson et al. 2020; 
Wielgus 2021; Broderick et al. 2022b). Motivated by this strong 
physical prediction, we proceed to fit a thin geometric cres- 
cent model with an asymmetric azimuthal brightness distribution 
combined with the standard splined raster and large-scale Gaus- 
sian model described above, collectively referred to as a "Hybrid 
Themage". A more detailed explanation of the model parameters 
and prior distributions is given in Appendix J. 

There have been important discussions about how well the 
Hybrid Themage model recovers the properties of the n = 1 pho- 
ton ring and by how much the recovered parameters are biased 
by contributions from the n = 0 ring (Tiede et al. 2022b). For 
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this reason, we cannot assert that this thin geometric component 
represents the n = 1 photon ring. Any model with an explicit ring 
component is likely to fare well when applied to M 87* EHT data 
with sufficient coverage. In a future paper, we intend to conduct 
a series of detailed synthetic data tests to better understand the 
behavior of the thin ring component. 

For the Hybrid Themage model, we fit the same scan- 
averaged, network-calibrated data we use in the standard 
Themage raster fitting, along with the same priors on the station 
gain parameters. We used the same sampler and adaptive tem- 
pering scheme to fit the Hybrid Themage models we also used to 
fit the standard Themage models. This type of image model has 
consistently produced good fits when applied to M 87* data, both 
in previous papers analyzing the 2017 data (Lockhart & Gralla 
2022; Broderick et al. 2022a) and in the current data set as seen 
in Appendix K. 

As mentioned above, we also repeat the same procedure done 
in M87* 2017 VI to directly model physical parameters in the 
visibility data. We used the THEMIS xs-ringauss model to pro- 
vide direct continuity with the 2017 EHT analysis, which was 
also the best performing crescent model featured in Sect. 6.1. A 
detailed description of the construction of the xs-ringauss model 
is described in Sect. 5.1 of M 87* 2017 VI, but we also briefly 
summarize the main components here. 

The primary ring was constructed by subtracting a uniform 
disk offset from a larger uniform disk, creating a “top-hat’ cres- 
cent, then "slashed" to create a linear brightness gradient. To 
quantify the depth of the central flux depression, we added a cir- 
cular Gaussian to act as a "floor" in the crescent center. We also 
added an additional asymmetric Gaussian pinned to the inner 
edge of the crescent at the point where the width is largest, 
inspired by the model from Benkevitch et al. (2016). The overall 
"image" was then smoothed with a Gaussian kernel, and two 
additional "nuisance" asymmetric Gaussian components were 
introduced to help model diffuse emission near the ring. A car- 
toon describing each model component is shown in Fig. N.1. 

The main difference between the Hybrid Themage and the 
xs-ringauss model construction is in the treatment of the dif- 
fuse emission. The xs-ringauss model used a pair of asymmetric 
Gaussians to help model any non-ring emission, and the Hybrid 
Themage used the adaptive raster components. The ring model 
in the Hybrid Themage was constructed using the same model 
object in the THEMIS code as the xs-ringauss, but with restricted 
priors to suppress the “pinned” Gaussian and enforce a small 
ring width. While the Hybrid Themage utilized a more flexi- 
ble method to model the non-ring emission, the xs-ringauss was 
more agnostic to the precise structure of the ring. 

We compare the xs-ringauss model to the scan-averaged 
network-calibrated data, similar to how we fit the raster Themage 
models. We utilized the same parallel tempered Hamiltonian 
MCMC sampler to explore the likelihood space for the xs- 
ringauss. We also fit the complex station gains in the same way as 
the raster models, using the same priors for the gain amplitudes 
and phases. We describe the detailed priors on the xs-ringauss 
parameters in Appendix N. 


6.2.2. Comrade direct modeling method 


We also define a fiducial model based on the best-performing m- 
ring model in Fig. 14 for analysis with Comrade. The fiducial 
model was taken to be an m-ring of order 2 that was allowed 
to have an asymmetric stretch. We included a flat emission 
floor in this model and constrained it to fill the interior of the 
m-ring. The resulting combination was convolved with a circu- 


lar Gaussian. This construction defines a nine-parameter, ellip- 
tical, floored m-ring model, which we refer to hereafter as an 
"mF-ring". An additional pair of elliptical nuisance Gaussians 
were also included to improve the fit quality. A summary of the 
mF-ring and its defining parameters is given in Fig. N.1, with 
additional details in Appendix N. We have also fit the fractional 
stretch of the mF-ring and the xs-ringauss compact flux density. 
We do not report a compact flux density for the mF-ring since 
we fit this model to closure quantities which are insensitive to 
the total flux. 


6.3. Direct modeling results 


Similar to the imaging methods in Sect. 5, each of these direct 
modeling methods fits the data well. We compare the model per- 
formance against the data in Fig. 15, where we show the com- 
plex visibility residuals for the Hybrid Themage and xs-ringauss, 
as well as the log closure amplitude and closure phase residu- 
als for the mF-ring. The gray points represent the scan-averaged 
data, the colored points represent the maximum likelihood model 
for the THEMIS-based models, and the maximum a posteriori 
model for the mF-ring. We find that each model performs very 
well when compared against the data, with small, unstructured 
residuals. Appendix N also presents image-domain representa- 
tions for each day and all observing bands, demonstrating con- 
sistency with the other imaging methods. 

Table 11 gives a summary of the physical model parameters, 
and Fig. 16 shows the model parameters and their comparison 
with the same model parameters from 2017. The most notewor- 
thy change from 2017 is that of the brightness asymmetry posi- 
tion angle (9 in Table 11), which has shifted from 160° to 212° 
as fitted by the mF-ring, and 202° as fitted by the xs-ringauss. 
The Hybrid Themage finds a position angle of about 215°. The 
slight discrepancy between the three models is likely due to our 
models definition of position angle, which is ambivalent to any 
flux contribution from the nuisance Gaussians and raster compo- 
nents. We note that image domain fits of the best-fit xs-ringauss 
and mF-ring models produce consistent position angle extrac- 
tion and thus conclude that the measured position angle from the 
mF-ring and the xs-ringauss model are consistent. 

We also measure a slight increase in diameter from 
43 uas as measured by the direct modeling methods in 2017 
(M 87* 2017 VI) to a median of 44.6 uas and 44.9 uas for the xs- 
ringauss and mF-ring models, respectively. The Hybrid Themage 
has a mean diameter even higher, at 45.5 uas, but the known 
inverse correlation between diameter and width for thin rings 
(M 87* 2017 VI) can partially explain this. It is expected that 
most of the observed 1.3mm emission originates from the 
inner portion of the accretion flow-in a region near the vicin- 
ity of the photon sphere (M 87* 2017 V; Dexter et al. 2012; 
Emami et al. 2023)-though this emission morphology produces 
images that are slightly larger than the predicted size of the black 
hole shadow. Moreover, this emission is subjected to vary in 
size resulting from accretion flow physics (Tiede et al. 2022b). 
Section 5.3 of M 87* 2017 VI gives an estimate of the theoreti- 
cal scaling and uncertainties associated with inferring physical 
features from the fitted parameters of the geometric models. 

The fractional stretch of the mF-ring allows for a mea- 
surement of image ellipticity which could originate from 
many possible features such as the shadow (Falcke et al. 2000; 
Cunningham & Bardeen 1973) or inner shadow (Chael et al. 
2021). However, the ability to extract this feature from 
data can be complicated by gaps in the EHT (u,v) cov- 
erage. It thus may not be directly related to any intrinsic 
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Fig. 15. Complex visibility (V) fit comparisons and normalized residuals for the Hybrid Themage and xs-ringauss, and the log-closure amplitude 
(In Ac) and closure phase (Wc) fit comparisons and normalized residuals for the mF-ring model. Filled points correspond to the real components 
of the complex visibility and open points correspond to the imaginary components. The colored points come from the maximum likelihood model 
for the Hybrid Themage and xs-ringauss, and from the maximum a posteriori model for the mF-ring. The gray points represent the scan-averaged 


data, with 1% fractional systematic noise added to all baselines. 


image geometry (Tiede et al. 2022c). The depth of the cen- 
tral brightness depression and the fraction width measure- 
ments in 2018 are consistent with their respective measurements 
from 2017. 


7. Discussion 


In this work we utilize the EHT observations from April 2018 
to obtain a new image of M87", building upon the findings 
from the 2017 EHT campaign (M87* 2017 I; M 87* 2017 IL; 
M 87* 201710; M87*2017IV; M87*2017 V; M87* 2017 VI; 
M87* 2017 VII; M 87* 2017 VIII). As shown in Fig. 17, repre- 
sentative images from band 3 April 21 data produced by both the 
imaging pipelines and direct modeling methods reveal the exis- 
tence of a ring-like emission structure around the central SMBH 
M 87*. This result across various independent methods demon- 
strates the robustness of our conclusions. The consistency of a 
ring-like emission structure on a 1-year time scale supports the 
interpretation of this structure as the black hole shadow of M 87* 
(M 87* 20171; M87* 2017 VI), consistent with the predictions 
from general relativity. 
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In this section, we discuss the main implications from the 
new images in 2018. We start by discussing the significance of 
the new diameter estimate and compare it to the previous diame- 
ter and mass estimate from the analysis of the 2017 data. We do 
not produce a new mass estimate, since this requires a dedicated 
calibration exercise using GRMHD simulations (M 87* 2017 VI; 
Sgr A* 2017IV; Sgr A* 2017 VD, which will be featured in a 
follow-up paper. Instead we use the scaling factor from the 
2017 analysis to construct a proxy for the gravitational diame- 
ter (M 87* 2017 VI), allowing us to compare the mass estimate 
from 2017 with the diameter measurements from 2018. Then 
we discuss the significance of the counter-clockwise shift in the 
position angle of the brightness asymmetry in the ring from 
2017 to 2018. Finally, we investigate the challenges related to 
estimating the compact flux density of M 87* on horizon scales 
with the sparse (u, v) coverage of the EHT. 


7.1. Consistency of the ring diameter 


The asymmetric ring's estimated characteristic parameters are 
provided in Tables 7-10 via image domain feature extraction 
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Fig. 16. Fitted features of the xs-ringauss, mF-ring, and Hybrid Themage fiducial models to the 2018 M 87* HOPS data. We include fits from 
bands 1 through 4 on April 21 and April 25. Each point shows the median value of the posterior distribution and the error bars indicate the 68% 
posterior probability range centered around the median. The blue line and band represent the median and 6896 confidence interval for the posterior 
generated by combining all bands and all days for the mF-ring model, while the red band is the equivalent for the xs-ringauss. The pink lines and 
bars represent the statistical mean and standard deviations of the Hybrid Themage method. The black line indicates the median over all days and 
bands of the 2017 M 87* analysis from xs-ringauss fits. The hashed region is the 6896 posterior probability interval taken from the 2017 M 87* 


xs-ringauss fits. 


and Table 11 via the direct modeling results. The parameters 
are broadly consistent across different calibration schemes and 
reconstruction methods. There is a tendency for the ring diam- 
eters extracted with the image domain feature extraction to be 
slightly smaller than those from the direct modeling methods. 
This tendency is also observed in the 2017 results and is related 
to the effective resolution of each method (M 87* 2017 VI). 

We combine the results from image domain feature extrac- 
tion and direct modeling to construct a median diameter across 
all methods (Tables 7-11). We first construct a diameter his- 
togram for each method, normalize the histograms to equally 
weigh each method, then combine the normalized histograms 


to produce a single diameter histogram across all imaging and 
direct modeling methods. We present the resulting median diam- 
eter and 68% intervals in the last column of Table 12. This lets us 
more easily compare the ring diameter of 2018 EHT campaign 
with that of 2017 (see Table 1 of M 87* 2017 I). We find the dif- 
ference in median diameter between these two campaigns is on 
the order of only ^1 uas, with error bars that almost completely 
overlap, indicating a persistent emission structure. 

Since the diameter estimates are model-dependent, we can 
instead convert the diameter to a common physical scale. For 
black holes, a natural quantity is the angular gravitational radius 
Og. We can construct this quantity by using the diameter values 
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Table 11. Summary of the direct modeling parameters. 


Data set Parameters 
Day Band Code Model d (uas) p d(deg) logis (fe) CF Jd 
April21 bl THEMIS  xs-ringauss  44.3+988 0.36'007 208.3730 = - 1.144928 1.017007 - 
Comrade mF-ring ang 0.13400% 210552. 1492 0.114993 
THEMIS Hybrid Themage — 46.602? 218.843! 1,114003 
b2 THEMIS xs-ringauss = 45.27 89 0.271903 206.2772 —1.41+936 1.027005 n 
Comrade mF-ring 4551919 04005 2100 50025 - 0,0209 
THEMIS Hybrid Themage | 46.515; 2164715 0.64*003 
b3 THEMIS xs-ingaus 43.622 0.227001 203.312  -2.12*04$ 1.03700 m 
Comrade mF-ring 44.5700 0.1904 212.6731. —1.817028 TT 0,0500. 
THEMIS Hybrid Themage — 44.9709, 216.177 Dau 
b4 THEMIS  xs-ringauss 444729 08:00  198.6*!5  —242*05  1.02*005 - 
Comrade mF-ring 44.6*017 Cts 212.1498 —1.47+9:% s 0.05+0:03 
THEMIS Hybrid Themage 45.2*?-02 210.743 05870 
April25 bli THEMIS xs-ringauss 43.81133 0.39'000 191.5731  —0.58:075 0.947007 = 
Comrade mF-ring 4662. O19, 207.91723  —0.81+920 obs 0.06+9-06 
THEMIS Hybrid Themage 46.0+9:88 si 224.71300 I 1.05*005 
b2 THEMIS xs-zingaus 43.7512 0.377008 1944759  —0.84:031 0.9700) E 
Comrade mF-ring 47.94733 0.14t007 2187 2". Klee iis 0.06*007 
THEMIS Hybrid Themage | 44.9*1 TT 2124153 is 1.09+0.05 
b3 THEMIS xs-ringauss 44.606 0.271903 207.7523 -1.351933 0.97006 i 
Comrade mF-ring Mor 029979 ITO", -d30705 T 0.08*002 
THEMIS Hybrid Themage — 45.502, T 214.3437 sles 1.06107 
b4 THEMIS  xs-ringauss 44.570533 0.23'006 199822  —122:931 0.987007 " 
Comrade mF-ring 45499 0234 210175 120922 0073094 
THEMIS Hybrid Themage | 44.703 206.9*22 LII 


Notes. d is the model diameter, f; is the model fractional width, ¢ is the model position angle, log mi fc) is the logarithm of the fractional central 
flux depression, CF is the model compact flux, and f is the mF-ring fractional stretch. We report the median values and the 68% confidence 


intervals. 


Table 12. Comparison of the extracted ring parameters for 2017 and 2018. 


2017 


2018 


Ring diameter 
Orientation position angle 


42 + 3 uas 
150° — 200° east of north 208° + 5° east of north 


43.3*15 uas 


(d) obtained via image reconstruction and direct modeling along 
with a scaling factor (o) related to GRMHD simulations: 


G) 


Since we know the true angular gravitational radius in the 
GRMHD simulations, we can produce image reconstructions of 
GRMBHD synthetic data and directly determine the scaling factor 
between the reconstructed diameter and the true 6,. This scaling 
factor is a function of the GRMHD black hole spin, inclination, 
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and temperature ratio for the electron-to-proton coupling in both 
the weakly and strongly magnetized regions. The scaling factor 
also depends on the particular image reconstruction model, and 
so in principle it should be unique for each method. 

In our previous 2017 analysis, we derived «œ values between 
10.7 and 11.5 with associated errors of ~10% across all methods, 
from imaging to direct modeling (M 87* 2017 I; M 87* 2017 VD. 
Because the œ value has some dependence on the observa- 
tional properties such as coverage and S/N, and because we 
use some methods in 2018 that were not used in 2017, it is 
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Fig. 17. Representative band 3 images of M 87* on April 21 from each imaging pipeline (top row). Fiducial images are shown for DIFMAP and 
the RML methods. The DIFMAP image is restored with a circular 20 was beam, as shown by the circle in the lower right corner. For THEMIS, 
Comrade, Hybrid Themage, the xs-ringauss, and the mF-ring a random sample drawn from the posterior is shown. Representative band 3 images 
of M 87* on April 21 from each imaging pipeline after blurring them with a circular Gaussian beam are shown on the bottom row — the FWHM 
of each beam is shown with the horizontal bar in the bottom right of each image. The eht-imaging, SMILI, THEMIS, Comrade, and Hybrid 
Themage images have been restored with 16.9, 16.1, 19.5, 18.8, and 14.2 uas FWHM Gaussian beams, respectively, to match the resolution of 
the DIFMAP reconstruction, whereas the xs-ringauss and mF-ring images were restored with a 20 uas FWHM Gaussian beam. The vertical dashed 
line separates the DIFMAP and RML methods from the Bayesian methods, and the solid vertical line separates the imaging methods which do not 
assume a ring-like structure, from the direct modeling methods which do assume a ring-like structure. 


necessary to re-calibrate the 2018 methods to determine a val- 
ues for a proper 2018 6, estimate. This would require image 
and model reconstructions of a large number of GRMHD syn- 
thetic data sets. This is a significant undertaking, and since 
this paper is primarily focused on presenting the first 2018 
images, it is beyond the scope of this work to produce a 2018 6, 
calibration. 

Nevertheless, it is still valuable to produce some compar- 
ison with the 2017 results, and to that end we construct a 
proxy for Og, d/a2917, which combines the diameters we mea- 
sure in 2018 with the œ values from 2017. For the 2018 meth- 
ods which have a 2017 6, calibration, we will use the corre- 
sponding 2017 a values from M 87* 2017 VI, either in Table 4 
for the xs-ringauss or Table 6 for DIFMAP, eht-imaging, and 
SMILI. For the 2018 methods which do not have a corre- 
sponding 2017 60, calibration, we will adopt œ = 11.0, which 
is the median value across all 2017 methods (see Table 1 in 
M 87* 2017 D). 

We show a comparison between the 2018 April 21 band 3 
d/a591; and 2017 6; values in Fig. 18 for each method. We show 
the 2018 points with colored diamonds and the 2017 points with 
black squares. Red (REx) and blue VIDA points correspond to the 
diameters estimated via image domain feature extraction, and 
the pink, purple, and dark blue points correspond to the Hybrid 
Themage, xs-ringauss, and mF-ring direct modeling methods. 
The uncertainty in the 2018 points is directly related to the 68% 
confidence intervals of the April 21 band 3 diameter estimates 
(see Tables 7 and 11). For the 2017 points, the uncertainty related 
to the 2017 observational systematics are slightly larger than the 
measurement uncertainty (see Sect. 7 of M 87* 2017 VI), and 
is represented by the solid black error bars. The uncertainty 
related to the diversity of GMRHD models used in the calibra- 
tion set dominates over all other uncertainties, and is shown by 
the dashed gray error bars. 

While the diameter uncertainties for some of the 2018 recon- 
structions are larger than the observational uncertainties in 2017, 
all methods are broadly consistent with each other, and all live 


within in the GRMHD uncertainty. This supports our previous 
interpretation that the ring-like structure in our images is consis- 
tent with the shadow of a supermassive black hole with a mass 
of ~6.5 x 10? Mo. Any discrepancies may disappear in a proper 
2018 6, calibration. 

We note that there are at least two new integral field spec- 
troscopic stellar kinematics mass estimates for M 87* recently 
published. One of these adapts new triaxial orbit models to 
Keck II telescope observations (Liepold et al. 2023) and favors a 
mass for M 87* of 5.3 x 10? Mo. The other estimate uses adap- 
tive optics observations from the VLT instruments MUSE and 
OASIS (Simon et al. 2024) and favors a mass of 8.7 x 10? Mo. 
However, when considering a different stellar mass profile in 
the inner region (and thus a different mass-to-light ratio), this 
method obtains a mass estimate of 5.5 x 10? Mọ. These ~5 x 
10? Mo mass estimates are within 1.5c of the 2017 EHT mass 
estimate. Arguably, there are significant systematic uncertainties 
in the non-horizon scale mass estimates of M 87* that require 
continued investigation. 


7.2. Position angle variation 


It is worth noting that there is a significant counter-clockwise 
change in the position angle of the brighter region in the ring- 
like structure from 2017 to 2018 (—30?). The shift direction is 
consistent with the prediction reported in M 87* 2017 V if the 
alignment of the large-scale jet is normal to the disk. 

The work by Wielgus et al. (2020) shows evidence of year- 
scale position angle variation of the ring's peak brightness, 
based on simple modeling of prototype EHT data. The present 
study is the first case where such a year-scale variation is 
unambiguously confirmed in the image domain (see Fig. 19). 
Furthermore, recent long-term monitoring studies of M 87* 
using longer-wavelength VLBI found a systematic position 
angle oscillation of the parsec-scale jet (Walker et al. 2018; 
Cui et al. 2023), which could be caused by flow instabilities 
(Walker et al. 2018) or precession of the central compact source 
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Fig. 18. Comparisons of d/a 9,7, which serves as a proxy for @,, using multiple methods for 2017 (squares) and 2018 April 21 band 3 (diamonds). 
For the 2018 methods which have a corresponding 6, calibration from the 2017 analysis (DIFMAP, eht-imaging, SMILI, xs-ringauss), we use the 
method-specific 2017 scaling factor to determine the 2018 d/a2017 values. For DIFMAP, eht- imaging, and SMILI we use the scaling factors from 
Table 6 of M 87* 2017 VI, and for the xs-ringauss we use the scaling factors from Table 4 of M 87* 2017 VI. For the 2018 methods which do not 
have a 2017 6, calibration (THEMIS, Comrade, mF-ring), we use a = 11.0, coming from the median a across all 2017 methods (M 87* 2017 I). 
The error bars for the 2018 points are representative of the 68% confidence intervals of the model-specific diameter estimates. The two image 
domain feature extraction methods are shown with red points for REx and blue points for VIDA. For the 2017 points, we show the measured 6, 
(black squares), the uncertainty due to differences in the 2017 observational details (74s, solid black error bars), and the uncertainty due to the 
diversity of the 2017 GRMHD library (my, gray dashed error bars). The gray horizontal line and shaded region represent the 2017 6, value and 
Otny uncertainty for the xs-ringauss model. While the uncertainty in diameter for some of the 2018 methods is larger than that for the same methods 
in 2017, the uncertainty related to the different GRMHD models used as the calibration set dominates and spans all methods. 
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Fig. 19. Comparison of the brightness position angle measured in the EHT observations during 2009—2018. The orange shadow covering (288 + 
10)? indicates the observational position angle range of the milliarcsecond-scale jet (Walker et al. 2018; Cui et al. 2023). The red points for 2018 
are the median values after combining the posteriors of all bands on April 21 and April 25 (see Table 11). The blue points for April 6 and April 
11 in 2017 are adopted from M 87* 2017 IV, M 87* 2017 VI and Wielgus et al. (2020, see the Table 4). The gray points for 2009-2013 are adopted 
from Table 3 in Wielgus et al. (2020), and represent the 68% confidence intervals. The posterior shapes for the 2009—2013 are non-Gaussian, and 
exhibit large tails. 


(Cui et al. 2023). M 87* exhibits year-scale morphological vari- and also the lack of inter-year EHT images. Further accumula- 
ations and a counter-clockwise shift from 2017 to 2018 at uas tion of EHT images over forthcoming years, along with parsec- 
scales. Nevertheless, whether the time variation we see in M87" scale jet monitoring, is required to understand better the origin of 
and its jet are physically linked with each other or not, there the year-scale variation of the ring-like structure and its possible 
could be some biases, given the large spatial gap between them connection to the large-scale jet. 
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Fig. 20. Distribution of best-fit position angle (in degrees) of the forward jet by fitting model images with the 2017 (April 6, high-band) and 
2018 (April 21, band 4) observations. The best-fit 10% of images (solid lines) among all (^18 000) images (dashed lines) are also shown. For 
reference, the vertical line shows the position angle ~288° of the milliarcsecond-scale jet (Walker et al. 2018; Cui et al. 2023), with the orange 
shadow area from (288 — 10)? to (288 + 10)°. While the fitted model images include different black hole spins, accretion types, and different 
electro-thermodynamics (M 87* 2017 V; M 87* 2017 VIII), the black hole spin vector is pointing away from Earth in all images. The 2018 EHT 
results are consistent with a black hole spin vector pointing away from Earth. 


In addition to the year-scale position angle variation of 
the brightest spot in the asymmetric ring (Wielgus et al. 2020), 
GRMHD simulations for the accretion environment around 
M 87* also show that the position angle of the brightest location 
of the asymmetric ring can vary due to the turbulent, magnetized 
accretion environment, with a time-scale much smaller than the 
observational cadence between 2017 and 2018 EHT observa- 
tions (see also Fig. 6 of M 87* 2017 V). As will be discussed 
further in a forthcoming paper, it is possible to apply the 2017 
and 2018 EHT observation results as independent constraints for 
GRMHD models of the black hole in M 87. 

Figure 20 presents the position angle distribution of the 
fitted jet directions, assuming the black hole spin vector is 
pointing away from Earth. To this end, from the image library 
applied in M 87* 2017 V and M 87* 2017 VIII, we select a group 
of model images which have the black hole spin axis point- 
ing away from Earth, then fit all the selected images to the 
2017 (April 6, high-band) and 2018 (April 21, band 4) EHT 
data by applying the snapshot model (SSM) scoring proce- 
dure introduced in M 87* 2017 V and M 87* 2017 VI. For a sum- 
mary of the model image library preparation in M 87* 2017 V 
and M 87* 2017 VIII), the one-fluid GRMHD simulations are 
initialized with a weakly magnetized torus of plasma orbit- 
ing in the equatorial plane of the black hole, and the rota- 
tional axis of the torus is aligned with the black hole axis. 
After the system evolves to a steady state, we perform radiative 
transfer to compute the thermal synchrotron emission to gen- 
erate images with the parameterized electron thermal dynam- 
ics (see M 87* 2017 V, for details). The selected models were 
prepared by the code iharm (Gammie et al. 2003), and the 
model images including different black hole spins?, accretion 
flow types, and assumptions for the ratio between the elec- 
tron and ion temperature in the simulation (see M 87* 2017 V, 
for details). 

Constrained by the 2018 EHT observation of M 87*, the dis- 
tribution of the forward jet direction in the model images is con- 
sistent with the observed milliarcsecond-scale jet (Walker et al. 


? The dimensionless spin parameter a, and inclination angle i between 
the observer's line of sight and the spin axis of the accretion flow are: 
a, = (—0.94, —0.5) for i = 17°, and a, = (0,0.5,0.94) for i = 163°. 


2018; Cui et al. 2023), which is ~288° counter-clockwise from 
the north. This supports the interpretation that the black hole 
spin vector points away from Earth, which was also strongly 
favored by the first M 87* EHT results (M 87* 2017 V). Addition- 
ally, with the model images considered, the position angle distri- 
bution of the forward jet direction in the 2018 observations has 
its peak closer to ~288° (vertical dashed line in Fig. 20). Assum- 
ing the large-scale jet is aligned perpendicular to the accretion 
disk (M 87* 2017 V), with many yearly observations to accumu- 
late a statistically significant number of independent images, we 
expect to see the position angle of the brightest region to be more 
similar to that seen in 2018 than in 2017. Since the timescale 
for position angle shifts in GRMHD simulations is small com- 
pared to the yearly cadence of EHT observations, we do not nec- 
essarily associate this counter-clockwise shift in the brightness 
position angle with a global rotation of the accretion flow, but 
instead consider the 2018 observations to be a snapshot of the 
most common orientation of the accretion flow. Additionally, we 
will also directly score a new GRMHD library against the 2018 
data using similar procedures to the ones found in M 87* 2017 V; 
M 87* 2017 VIII; Sgr A* 2017 V, and leverage the 2017 and 
2018 data to score the GRMHD simulations across multiple 
years. 


7.3. Compact flux density estimates from Bayesian imaging 
methods 


The EHT has very few baselines that probe the intermediate 
region in (u, v) space between the intra-sites and about 2 GA. The 
only baseline that probe this region of (u, v) space is LMT-SMT. 
Even though the LMT has a very sensitive 50 m dish, it typi- 
cally experiences relatively large pointing errors, increasing the 
recovered station gain amplitudes. While the LMT seems to have 
much better pointing in 2018 than in 2017, we have about half 
the scans in 2018 compared to 2017 and no LMT participation 
on April 25. This limited data means it is much more challeng- 
ing to constrain the compact flux density in M 87* during the 
2018 EHT observation campaign. The longer-wavelength and 
zero-baseline analysis described in Sect. 4 are consistent with 
the findings from the Bayesian methods: there is a large amount 
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Fig. 21. Model compact flux density estimates from the Bayesian Image 
reconstruction methods for the April 21 HOPS data. The error bars rep- 
resent the 68% confidence intervals around the median (square points). 


of uncertainty in what we should expect the compact flux density 
to be at EHT-scales. 

To illustrate this, we investigate the compact flux density 
estimates produced by the various Bayesian methods that fit 
the complex visibilities (THEMIS) or the visibility amplitudes 
(Comrade). For the THEMIS raster and Hybrid Themage mod- 
els, we quantify the compact flux density by summing the flux 
density from the raster and the raster plus crescent, respectively. 
For Comrade, the compact flux density is a model parameter 
which sets the total flux density of the raster as described in 
Sect. 5.1.5. The THEMIS xs-ringauss model constructs the com- 
pact flux density estimate by summing the flux density from the 
crescent and nuisance Gaussians. 

We summarize the compact flux density estimates across all 
observing bands on April 21 from the THEMIS raster, Hybrid 
Themage, Comrade raster, and THEMIS xs-ringauss models in 
Fig. 21, using the HOPS data. There are two clear modes: a “low- 
estimate" for compact flux density around 0.6 Jy, and a “high- 
estimate” for compact flux density around 1.0 Jy. The THEMIS 
raster model is the only model that consistently produces com- 
pact flux density estimates below 0.7 Jy for all bands. In com- 
parison, the THEMIS xs-ringauss and Comrade raster produce 
compact flux density estimates that are greater than 0.9 Jy. The 
Comrade raster prefers larger compact flux density because of 
the larger FOV and the choice of prior distribution. We found 
that the compact flux density strongly depends on the Comrade 
prior assumptions because of the lack of intermediate baselines. 
The Hybrid Themage model prefers lower compact flux densities 
in band 3 and band 4, but a higher compact flux density in band 
1, with uncertainty that spans 0.5 Jy. The error bars in Fig. 21 are 
taken from the 68% confidence intervals and that, generally, the 
methods that produce low compact flux density estimates do not 
overlap with the high compact flux density estimates. 

There is no consensus among the Bayesian methods about 
what the compact flux density should be in the 2018 EHT data. 
The chain parameters for the Hybrid Themage reconstruction 
of the April 21 band 1 data give us clues about how we might 
achieve such high values for compact flux density. In Fig. 22, 
we show the total compact flux density, the raster FOV, and the 
crescent flux density. Recall that for the Hybrid Themage, we 
construct the compact flux density by adding the crescent flux 
density and the flux density from all the raster points. We see 
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Fig. 22. Triangle plot comparing the compact flux density (CF), raster 
FOV in the x and y directions, and crescent flux (Jy) parameters for 
the Hybrid Themage reconstruction of April 21 band 1. The compact 
flux density and raster FOV parameters exhibit clearly bi-modal distri- 
butions in this parameter space. The color regions represent the 99%, 
90%, and 50% quantile regions. 


clear bi-modal distributions in the compact flux density, FOV,, 
and to a lesser degree, in the crescent flux density. This sug- 
gests that, at least for this reconstruction, adding the entire raster 
flux density may not be appropriate for computing the compact 
flux density; the raster is elongated in one direction and places 
some flux outside of what we might consider the compact region. 
This is evident in the April 21 band 1 Hybrid Themage image in 
Fig. N.2, where the raster extends flux in the northwest direction. 
The most crucial parameters of interest, such as the diameter 
and the position angle, do not seem to depend strongly on the 
recovered compact flux density. This is consistent with the find- 
ings from the RML and CLEAN methods, which are not able to 
independently constrain the ratio between the total flux density 
and the compact flux density. The measured diameters and posi- 
tion angles are entirely consistent between models that produce 
different values of compact flux density (see Appendix M). 


8. Conclusion 


We present the results of the 2018 EHT observations of M 87* 
at 1.3 mm in order to further investigate the nature of the black 
hole shadow, and the year to year variation in the image struc- 
ture. During the 2018 observations, the GLT participated for 
the first time as part of the VLBI array and provided addi- 
tional baselines to improve the (u,v) coverage. Following the 
example of the 2017 data analysis, we use multiple independent 
calibration, imaging, and modeling methods to analyze and 
interpret the 2018 EHT data. We analyze data on two indepen- 
dent epochs (April 21 and April 25) across four independent 
frequency bands. The 2018 data show a clear null in the visi- 
bility amplitudes at 3.4GA, and exhibit significant variation in 
the time-dependent closure phases compared to 2017. The dif- 
ferences in the closure phases indicate a significant change in 
the asymmetric structure from 2017 to 2018. 
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All of the imaging and modeling methods indicate the pres- 
ence of a ring-like structure, and the extracted ring characteris- 
tics are consistent across all bands and days. Similar to the 2017 
data, the current EHT data does not provide strong constraints 
on extended or diffuse emission outside the ring. All methods on 
April 21 data suggest a ring diameter of ~42 uas for the image 
domain analyses (Sect. 5.3) and ~45 uas for the direct modeling 
methods (Sect. 6.3). When combining all methods with equal 
weights, the median diameter of the ring in 2018 is 43377 uas, 
in good agreement with the diameter measured in the analysis of 
the 2017 data. 

The lack of variability in the ring diameter between 2017 
and 2018 is consistent with the predictions from general relativ- 
ity for strongly lensed emission around a black hole. In contrast 
with the 2017 data, we detect a significant shift in the position 
angle of the ring brightness asymmetry from ~180° as measured 
in 2017 to ~210° in the 2018 data. This shift in position angle 
is consistent with the expected variability from GRMHD simu- 
lations. In particular, when converting the brightness asymmetry 
position angle to a nominal black hole spin direction, the 2018 
image is more consistent with the orientation of the large-scale 
jet seen at longer wavelengths. 
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Appendix A: Flux density calibration parameters 
and uncertainties 


The visibility amplitudes need to be converted from correlation 
coefficients to units of flux density, considering the different tele- 
scope sensitivities across the array. This is achieved by multi- 
plying the normalized visibility amplitudes with the geometric 
mean of the derived system equivalent flux densities (SEFDs) of 
the two stations of each respective baseline. The SEFD is given 
by (e.g., M 87* 2017 IID: 

dipb m. 

Nph DPFU X gg 

where T3, is the effective system noise temperature, corrected 
for atmospheric attenuation. The degrees per flux density unit, 
DPFU, is a conversion factor from units of temperature (K) to 
units of flux density (Jy), and gz is the normalized elevation- 
dependent gain curve of the telescope. npn is the time-dependent 
phasing efficiency for a phased array, and has a value of unity 
for single dish telescopes. The DPFU can be estimated from the 
antenna aperture efficiency (ņa) as: 


(A.1) 


A eom 
DpFU = “= 


E (A.2) 


where Ageom is the geometrical area of the antenna dish in units 
of m?, and k = 1.38 x 10? is the Boltzmann constant in units of 
Jy m?/K. 

A summary of the station DPFUs and gain curve parame- 
ters for the EHT array in 2018 is provided in Table A.1, together 
with their estimated uncertainties. Details on their derivation for 
each station are described by Koay et al. (2023a). We expect 
the DPFU uncertainties to be the dominant source of systematic 
errors in the flux density calibration of EHT data. Ts values are 
not expected to vary much within the typical scan duration of a 
few minutes, so their uncertainties are expected to be no larger 
than a few percent. The GLT is an exception, due to sub-optimal 
performance in 2018 when it was still only partially commis- 
sioned, where T, uncertainties are also significant (Koay et al. 
2023b). 

We note in particular that the quoted errors for the GLT and 
LMT in Table A.1 are very likely to be underestimated, due to 
insufficient a priori data collected by the stations, as described 
in Appendix D. We therefore obtain an independent set of con- 
straints on the amplitude calibration uncertainties for these two 
stations in Appendix E. 


Table A.1. EHT station flux density calibration parameters and their uncertainties for the 2018 observing campaign. 


Station 


DPFU (K/Jy) gE 

LSB-RCP LSB-LCP USB-RCP USB-LCP uncertainties B Eo 
ALMA? 1.19 1.19 1.19 1.19 + 10% 0 0 
APEX? 0.0253 0.0262 0.0259 0.0270 + 5.5% 0.00002 + 3.6% 36.6 + 1.0% 
GLT* - - 0.00885 0.00885 + 7% 0 0 
LMT 0.188 0.126 0.188 0.126 + 22% 0 0 
SMT 0.0188 0.0188 0.0182 0.0182 + 6% 0.000082 + 10.4% — 57.6 + 2.096 
JCMT? 0.0296 - 0.0296 +(11-14)% 0 0 
IRAM-30m 0.144 0.138 0.142 0.137 x 1996 0.00018 + 5.3% — 43.7 € 1.396 
SMA* 0.040 0.040 0.040 0.040 +(5 - 15)% 0 0 
SPT 0.00698 0.00731 0.00698 0.00731 + 10% 0 0 


Note. DPFU values are quoted for both lower (LSB) and upper (USB) sidebands, and both polarizations (RCP and LCP). For phased arrays 
(ALMA and SMA), the DPFUs represent the combined sensitivity of all phased dishes. The gain curve parameters B and Ep as a function of 
elevation E are given based on a gg = 1 — B(E — Eo)? parameterization (see Janssen et al. 2019b; Koay et al. 2023a, for further details). "The 
ALMA DPFU uncertainty is based on the overall 10% systematic uncertainties associated with the QA2 flux calibration tables. "APEX measures 
DPFUS for each of the four bands separately. LSB and USB DPFUs shown here are the mean of the bands 1 and 2 DPFUs, and bands 3 and 4 
DPFUs, respectively. “GLT DPFU uncertainties are based on the scatter of three measurements, inclusive of flux density model errors, but may 


be larger due to significant antenna astigmatism. There is also a ~15% T$, 


, uncertainty contributing to the total SEFD uncertainty. For the JCMT 


DPFU, the nighttime value is shown. *The SMA DPFU uncertainty is based on the dominant 5—15^46 uncertainty on the phasing efficiency. 
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Appendix B: Non-closing systematic error budget 


The closure phases (Wc), log closure amplitudes (In Ac), clo- 
sure trace phases and closure trace log amplitudes 7, the lat- 
ter two of which are novel diagnostics described in detail by 
Broderick & Pesce (2020), are independent of station-based gain 
errors, but are sensitive to systematic non-closing systematic 
errors. Figure B.1 shows the normalized distributions of the dif- 
ferences in the closure quantities between bands (band 1 — band 
2, band 3 — band 4), polarizations (RR-LL), as well as the dis- 
tributions of the trivial closure quantities derived from triangles 
or quadrangles with at least two co-located stations. 

In an ideal scenario where polarimetric leakages are absent, 
the differences in closure quantities between bands and polar- 
izations, as well as the trivial closure quantities, is expected to 
follow a normal distribution with a mean value of zero. The pres- 
ence of non-closing systematic errors can lead to deviations from 
the expected normal distribution, as can be seen for the blue his- 
tograms in Fig. B.1. 

The non-closing systematic error, s, can thus be determined 
by solving for its value while enforcing the condition of: 


X 
mad, (<=) = mado ————|-21, 
o 


|) 2 
guts 


where œ is the total uncertainty associated with the random vari- 
able X, and c is the known a priori thermal uncertainty. In this 
construction, X is a random variable representing the closure 
products being examined, for instance band 1 — band 2 closure 
phase difference, RR-LL log closure amplitude difference, etc. 
mado is the modified median absolute deviation given by: 


(B.1) 


mado(Y) = 1.4826 med(]|Y]), (B.2) 
where med denotes the median, and the factor 1.4826 scales the 
result so that it acts as a robust estimator of standard deviation 
for a normally distributed random variable Y with zero mean. 

The non-closing systematic errors estimated for M87* and 
3C 279 using the various tests are presented in Table B.1. To 
avoid biases arising from low S/N, only closure quantities above 
the S/N > 6 threshold are used. 


Table B.1. Non-closing systematic uncertainties, s (and in units of thermal noise, s/o), for M87* and its calibrator 3C 279 estimated using various 


statistical tests on both CASA and HOPS products. 


Source Test CASA HOPS 
s s/o n s s/aq n 

M87* RR -LL closure phases 0.0? 0.0 339 00 0.0 446 
band 1 — band 2 closure phases 2.6? 0.3 229 0.0° 0.0 300 
band 3 — band 4 closure phases 1.5? 0.2 192 0.0° 0.0 251 
trivial closure phases 0.0° 0.0 2436 13° 0.2 500 
RR — LL log closure amplitudes 13.2% 0.7 158 40% 0.2 294 
band 1 — band 2 log closure amplitudes 5.4% 0.3 139 2.2% 0.1 255 
band 3 — band 4 log closure amplitudes 0.0% 0.0 116 00% 00 171 
trivial log closure amplitudes 0.0% 0.0 60 00% 0.0 119 
band 1 — band 2 closure trace phases 1.2? 0.2 43  0.0* 0.0 103 
band 3 — band 4 closure trace phases 3.4? 0.5 60 3.8 0.5 114 
trivial closure trace phases 0.0? 00 118 35 0.5 152 
band 1 — band 2 log closure trace amplitudes 0.0% 0.0 43 24% 0.2 103 
band 3 — band 4 log closure trace amplitudes 0.0% 0.0 60 53% 04 114 
trivial log closure trace amplitudes 4.8% 04 118 00% 0.0 152 

3C 279  RR-LL closure phases 1.8° 0.5 680 1.9° 0.6 687 
band 1 — band 2 closure phases 1.7? 0.5 404 L8 0.6 407 
band 3 — band 4 closure phases 1.9? 0.5 2416 L8 0.5 420 
trivial closure phases 0.9? 0.6 503 0.8° 0.6 503 
RR — LL log closure amplitudes 5.4% 0.7 833 54% 0.7 846 
band 1 — band 2 log closure amplitudes 5.6% 0.7 604 65% 0.9 605 
band 3 — band 4 log closure amplitudes 5.1% 0.7 619 64% 0.8 638 
trivial log closure amplitudes 3.2% 0.7 487 29% 07 491 
band 1 — band 2 closure trace phases 3.0? 0.9 378 38 1.2 393 
band 3 — band 4 closure trace phases 3.4° 10 359 3.» 0.9 412 
trivial closure trace phases 0.5° 0.3 316 0.5° 0.3 320 
band 1 — band 2 log closure trace amplitudes 8.5% 14 378 90% 16 393 
band 3 — band 4 log closure trace amplitudes 7.0% 1.1 359 6.3% 1.0 412 
trivial log closure trace amplitudes 2.7% 10 316 24% 0.8 320 


Note. Only data from Apr 21 and Apr 25 have been used. n is the number of closure quantities used for each test. 
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CASA : M87*, 3C279 
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Fig. B.1. Diagnostic plots showing the normalized distributions of various closure quantities of M87* and 3C 279 data from both the CASA (top 
four rows) and HOPS (bottom four rows) reduction pipelines. For each block corresponding to each pipeline, the first two rows show band 1 — band 
2 (first column), band 3 — band 4 (second column), RR-LL (third column), and trivial (fourth column) closure quantities. The bottom two rows 
show band 1 — band 2 (first column), band 3 — band 4 (second column), and trivial (third column) closure trace quantities. Only data from April 
2] and 25 have been used. The distributions prior to (blue) and after (red) accounting for the estimated systematic uncertainties, s, are shown. The 
values of s for each source and reduction pipeline are given in Table B.1. In the top left corner of each distribution, the number of >30 outliers 
are given considering thermal noise only, followed by the number of outliers considering thermal plus systematic noise for o in parenthesis. These 
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Appendix C: Temporal coherence 
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Fig. C.1. Cumulative histograms of the amplitude ratios between coher- 
ent averaging for entire scans (Agcan), and coherent averaging for 2s 
before incoherent averaging over scans (A»,), for both HOPS (blue) and 
CASA (orange) M87* data on both April 21 and April 25. The gray his- 
togram shows the corresponding ratios from the HOPS data on the same 
days with no atmospheric phase corrections applied. For each pipeline, 
the fraction of data with coherence above 90% (i.e., with decoherence 
losses of no more than 10%) when averaged over the full scan is indi- 
cated. 


We quantify the level of decoherence loss of the data calibrated 
using the HOPS and CASA pipelines using the ratio Agcan/A2s 
for each scan (M 87* 2017 IIT). Ascan is the debiased amplitude of 
the coherently averaged visibilities over the entire scan, and A», 
is the amplitude derived from data coherently averaged over 2s, 
then incoherently averaged over the entire scan. 

Figure C.1 shows the cumulative distribution of the quan- 
tity Agcan/A2s for 2342 unique scans (per band, polarization, and 
baseline) of the M87* April 21 and 25 data with S /N > 7, com- 
mon to both the HOPS and CASA pipelines. We find that 92.996 
of HOPS and 85.696 of CASA scan-averaged data show deco- 
herence losses of less than 10%. When selecting April 21 data 
only, 97.4% and 87.9% of HOPS and CASA scan-averaged data 
show better than 90% coherence. For April 25, 86.3% and 82.3% 
of HOPS and CASA scan-averaged data show better than 90% 
coherence. 

The level of coherence in the 2018 data is slightly lower than 
that of 2017 M 87* 2017 III, due to (i) overall poorer weather 
at some sites, including Chile, Hawai‘i, and at the SMT (espe- 
cially on April 25), as well as (ii) the poor sensitivity levels 
of the GLT, as described in Appendix E. Overall, CASA data 
have lower coherence levels than HOPS data, mostly for base- 
lines with lower S/N such as those including the GLT. For exam- 
ple, when we exclude GLT baselines, the percentage of CASA 
scan-averaged data showing decoherence losses of less than 1096 
increases by a few percent, and becomes comparable to that of 
HOPS. 


Appendix D: Data issues 


In this appendix, we provide details on data issues that were 
identified in the 2018 L1 data package pertaining to M87* (and 
3C 279), and steps taken to address these issues. Our analysis 


indicated that the flux densities on APEX baselines are under- 
estimated after a priori amplitude calibration, when compared 
to the flux densities of ALMA baselines to other stations. We 
also found 25% differences between RCP and LCP amplitudes 
on all APEX baselines in band 1, and a smaller but still sig- 
nificant 13% difference between the RCP and LCP amplitudes 
in band 4. Follow-up diagnostics found this issue to be caused 
by low power levels of the intermediate frequency (IF) signal 
and too large block-down converter attenuation settings for the 
ROACH2 Digital Back End (R2DBE; Vertatschitsch et al. 2015) 
associated with RCP in bands 1 and 4. These issues are corrected 
at the network calibration stage. There are also large amplitude 
drops on APEX baselines after UTC 01:00 on April 25, found to 
be caused by 180? peak-to-peak phase modulations whose origin 
is not yet determined to date. We therefore flagged and excised 
all the affected APEX data on April 25. The APEX baselines 
also show large amplitude jumps in several scans, typically at 
the beginning and end of a scan, which we also excised during 
post-processing. 

While the above problematic L1 data have either been cor- 
rected or excised during post-processing, there remain a number 
of issues in the data used in our analyses. The GLT was still 
under commissioning and participated on a best effort basis in 
2018. Due to astigmatism and a low ~ 22% aperture efficiency, 
the visibilities on GLT baselines have lower S/N compared to 
that of the other stations. The telescope sensitivity was also not 
well characterized in 2018, so large uncertainties in the ampli- 
tude calibration are expected (Koay et al. 2023b). Due to the low 
telescope sensitivity, pointing sources could not be tracked down 
to elevations of < 10°, resulting in a poorly constrained pointing 
model at those low elevations and additional amplitude losses 
for sources observed at those elevations. This does not signifi- 
cantly affect M87* which was observed at higher elevations, but 
severely compromises the quality of GLT data for 3C 279 which 
was observed between 6° to 8°. 

The abrupt halt of the LMT observations in the middle of 
the 2018 campaign meant that there were no antenna sensitiv- 
ity measurements recorded at the time. The amplitude calibra- 
tion applied uses sensitivity measurements obtained in 2020, 
based on the assumption that the antenna characteristics have not 
changed in the intervening period. Also there were focusing and 
pointing issues being corrected on the fly while observing M87* 
and 3C 279 on April 21, possibly leading to larger amplitude cal- 
ibration errors. No elevation dependent gain curves were applied 
owing to a lack of aperture efficiency measurements at low ele- 
vations, where large gain errors of 20% (or possibly higher) are 
expected. This is not expected to significantly affect the M87* 
and 3C 279 observations, since these targets were observed at 
high elevations of 50° to 77°. 

The SMA phasing efficiencies on April 25 are poor (mostly 
below 50% between UTC 03:00 and 05:30) and have a large 
scatter, resulting in a large scatter in amplitudes on SMA base- 
lines. The PV station also experienced poor weather on April 
25 after UTC 02:00, leading to a large scatter in amplitudes and 
phases on associated baselines. There are large drops in ampli- 
tudes on all ALMA baselines on April 21 between UTC 04:12 to 
04:17, attributed to phasing issues on ALMA. These are partly 
corrected during a priori amplitude calibration and network cal- 
ibration, but a residual 40% scatter or drop in amplitude remains 
on the corresponding scan, which can be fixed during model- 
dependent self-calibration downstream. 
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Appendix E: Additional LMT and GLT gain 
constraints from overlapping visibilities 


As mentioned in Appendix D, the GLT and LMT station sen- 
sitivities are not well constrained (see also Appendix A), and 
their residual gain uncertainties are expected to be larger than 
the amplitude calibration uncertainties provided in Table A.1 
in Appendix A. To obtain additional a priori constraints on the 
amplitude calibration uncertainties for baselines to these two sta- 
tions, we examine and compare amplitudes on the baselines that 
overlap with that of the LMT and GLT in (u, v) space. 

The ALMA-LMT and ALMA-SPT baselines overlap for the 
calibrator source 3C 279, (Fig. E.1 left panel). Assuming that 
the SPT sensitivity is well-characterized and that the amplitude 
calibration uncertainties are relatively small in comparison, we 
find the amplitudes on LMT baselines to be under-calibrated by 
~ 34% in bands 1 and 2, ~ 23% in band 3 (Fig. E.1, right panel), 
and ~ 16% in band 4. 


For the GLT, there are no overlapping baselines for M87*. 
Amplitudes on GLT baselines for 3C 279 cannot be used to con- 
strain its amplitude uncertainties for M87* due to the large point- 
ing offsets (see Appendix D). The closest baselines in (u, v) space 
where the amplitudes can be compared for M87* are between 
ALMA-LMT and GLT-PV, where we find the amplitudes on the 
GLT-PV baselines are lower than ALMA-LMT amplitudes by 
up to ~ 50% in bands 3 and 4. Considering that the LMT itself 
is under-calibrated, the amplitude errors for GLT are likely to be 
larger than 50%. 

These a priori constraints are very crude, zeroth order esti- 
mates, but provide useful consistency checks with downstream 
station-based gains derived from model-dependent calibration, 
that is self-calibration during the imaging analysis. They are also 
important to better characterize the station uncertainties used in 
the generation of the synthetic data in Appendix G, so that they 
are reasonably consistent with that of the actual M87* data. 
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Fig. E.1. Left: (u, v) coverage for 3C 279 observations in band 3 on 2018 April 21 (colored points), showing the overlap between ALMA-LMT 
and ALMA-SPT baselines (highlighted with a box). Right: band 3 flux densities on ALMA-LMT (blue) and ALMA-SPT (orange) as a function 
of the projected baseline length in the east-west direction, in units of wavelength, demonstrating the under-calibration of amplitudes on LMT 
baselines. HOPS data are shown in this figure, but are consistent with that from CASA. 
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Appendix F: Derivation of constraints on the total 
compact flux density and source size 


We estimate constraints on the total compact flux density by 
using the EHT data themselves, and by comparing the flux 
density from quasi-simultaneous multi-wavelength observations. 
These two methods are similar to those used in M 87* 2017IV 
and Sgr A* 2017 II. We provide the details below. 


F1. Constraints from EHT observations 


Here, we derive constraints on the total compact flux density and 
the source size. We use the procedure described in Appendix B.1 
of M 87* 2017 IV, which was updated in Sgr A* 2017 II to incor- 
porate our uncertainties for the antenna gains using the DPFU 
uncertainties. This procedure included four constraints for the 
previous M87* analysis. However, the 2018 data do not contain 
crossing points in the (u, v) space for the SMT-PV and LMT-PV 
baselines; therefore we cannot apply the 4th constraint. More- 
over, only the observations on April 21 and 22 contain the short 
intra-site baseline, SMT-LMT, which is necessary for this anal- 
ysis. Thus we can derive the three constraints summarized below 
only for the first two observing days. For reference, the total flux 
density Fi. estimated from the ALMA interferometric data is 
1.13 Jy for both April 21 and 22. This value is an average of the 
total flux values derived for each of the bands 1-4. 

Constraint 1. The first constraint is based on the fact that the 
visibility amplitudes on short baselines can be approximated by 
a circular Gaussian visibility function, 


enu") 


41n2 vey 


Vg(u; Vo, 0) = Vo exp [- 
where Vo is the total flux density of the Gaussian source, |u| is 
the length of the baseline, and 0 is its FWHM in radians. The 
intermediate-to-long baselines tend to measure larger correlated 
flux density than what is expected for an equivalent Gaussian 
source. 

We can deduce that the measured amplitude ratio of the 
ALMA-LMT over SMT-LMT baselines will be larger than the 
corresponding ratio from a circular Gaussian source model. Con- 
sequently, the FWHM size of a circular Gaussian determined by 
the amplitude ratio between SMT-LMT and ALMA-LMT base- 
lines provides an estimate of the minimum compact source size 
cpet that is not significantly affected by the intrinsic fine-scale 
source structure: 
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Here, *V;..; denotes the true visibility on the baseline i — j. Con- 
sidering the gain uncertainties, the ratio of the true visibility 
amplitudes is lower-bounded by 
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where V,- ; and Ag; denote the calibrated measured visibility and 
the gain uncertainty, respectively. We take the median of the vis- 
ibility amplitude ratio from the collection of constraints derived 
for each single VLBI scan to derive a robust estimate of the min- 
imum source size. 

Constraint 2. The second constraint comes from the curva- 
ture of visibility amplitudes between the intra-site baselines and 


(F.3) 


the LMT-SMT baseline. Since the compact flux density should 
not exceed the total flux density measured with the intra-site 
baselines, the amplitude drop from the intra-site to LMT-SMT 
baselines gives the maximum limit of that from the compact flux 
density to LMT-SMT baseline. Therefore, it gives the maximum 
limit of the source FWHM size with an equivalent circular Gaus- 
sian as 
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We assumed the uncertainty of the total flux Agi = 0.1. Simi- 
larly to Constraint 1, the median of the ratio is adopted to miti- 
gate the effects of statistical errors. 

Constraint 3. The minimum compact flux density can be 
derived by the maximum amplitudes of LMT-SMT baseline, 
since the visibility amplitudes are maximum at zero baseline 
length and the a priori calibration should not underestimate the 
station sensitivity. Therefore the compact flux density is con- 
strained as 


Fepet 2 (1 — Agi ner + TA |VtMmT-smTl. 


The equivalent circular Gaussian gives a stronger constraint with 
the minimum source size (Constraint 1) extrapolated from the 
LMT-SMT baselines, described by, 
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These constraints are validated using the various synthetic 
models described in Appendix G. The upper and lower limits on 
the source size obtained from Constraints 1 and 2 are along the 
directions of LMT-SMT and LMT-ALMA baselines. 

In summary, Constraints 1—3 lead to the conclusion that 
0.30Jy € Fa < 1.13 Jy on April 21 and 22. If we use Eq. F.7, 
the minimum compact flux is 0.38 Jy. The compact source size 
has an equivalent Gaussian FWHM ranging between 39 and 
98 uas on the first two observing days. 


F.2. Constraints from quasi-simultaneous multi-wavelength 
observations 


Supplementary, we also make use of quasi-simultaneous VLBI 
data at longer wavelengths obtained through part of a large 
MWL observing campaign coordinated with the EHT 2018 
observations. Our approach assumes that the total compact emis- 
sion of M87* at v « 230 GHz (A = 1.3mm) Fiot,230 originates 
from the combination of (i) a compact emitter with unknown 
flux density Fepct,230 and spectral index at 230 GHz, and (ii) a 
large and diffuse (milliarcsecond-scale) jet, which can be char- 
acterized by its total flux density Fjet230 = Fot230 — Fepet,230 
and a steep spectral index of aj4 ~ —(0.7 — 1.0), as measured 
at cm- and mm-wavelengths (with typical errors of ~ 0.3; see 
Hovatta et al. 2014; Hada et al. 2016). Then, by measuring Fjet 
and Fot at lower frequencies, and extrapolating them to 230 GHz 
using a power-law model, we can estimate the total compact flux 
density Fepet at 230 GHz: Fepet.230 = Ftot,230 — Fjet,230- 

We utilize the data obtained with the East Asian VLBI 
Network (EAVN; Anetal. 2018; Cuietal. 2021; Cho et al. 
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2022) at 22 and 43 GHz, the Korean VLBI Network (KVN; 
Lee et al. 2016) at 43, 86, and 129 GHz, and the Global Mil- 
limeter VLBI Array (GMVA)+ALMA at 86 GHz. All of these 
observations were conducted between 2018 March 25 and April 
21. For details of data reduction and calibration, see Cho et al. 
(2017), Cui et al. (2021) for EAVN, Lee et al. (2016), Kim et al. 
(2018b) for KVN, and Lu etal. (2023) for GMVA«* ALMA, 
respectively. 

The total flux densities of the nuclear region (i.e., Fi = 
Fiet + Fa) at 22, 43, and 86 GHz are measured from the EAVN 
and KVN images by integrating the flux densities over a large 
window of 500 x 500 jas? centered on the radio core. We obtain 
total flux densities of Fig»? > l.lJy, Fis43 œ% 0.9-1.1 Jy, 
Frogs & 0.9Jy and Frioti29 + 0.6-0.8 Jy at 22, 43, 86, and 
129 GHz, respectively. The KVN 129 GHz data suffer from 
larger uncertainties due to a stronger influence from weather. 
We fit a power-law model to the total flux densities (i.e., Fig; « 
y**w). We find aig ~ —0.17 when 129 GHz data are included, 
while a; ~ —0.11 when they are excluded. By extrapolating a 
power-law model with a range aig; ~ —(0.11 —0.17) to 230 GHz, 
we estimate Fiot,239 = 0.72 — 0.80 Jy. While this range of values 
is lower than the ALMA interferometric measurements during 
the EHT observations (~1.11—1.18Jy; Table 3 in Goddi et al. 
2021), it only includes the total flux density within the inner 
500 x 500 yas”, which is nearly three orders of magnitude below 
the resolution limit of ALMA. 

The compact flux density Fepct,s6 at 86 GHz is estimated 
from a GMVA+ALMA image that has a higher angular resolu- 
tion than the KVN images (Lu et al. 2023). We obtain Fonct.g6 ~ 
0.51 Jy by integrating over a smaller window of 100 x 100 uas? 
centered on the radio core. We can then derive a jet flux density 
at 86 GHz to be Fjet,86 = Fiot,86 = Fepct,86 z 0.39 Jy. 

Finally, we compute the expected Fepct,230 by estimating 
Fjet230 and subtracting it from the previously estimated F'ot,230- 
As explained above, we assume that Fjet,230 follows a power law 
that can be extrapolated up to 230 GHz (i.e., Fjet,230 ec v****). For 
a range of ag4 = —(0.7 — 1.0), we find Fje,230 = 0.15—0.20J y. 
This value gives a range of Fepct230 © 0.5 — 0.7 Jy within the 
central 100 x 100 uas? region. We emphasize that this estimate is 
based on lower-frequency data and assumptions on the jet spec- 
tral index. Therefore it should not be regarded as a tight con- 
straint on Fepct,230 but more as a useful reference. 


Appendix G: Imaging validation 


Here we describe the synthetic data generation process and sub- 
sequent reconstructions, as well as the resulting distribution of 
Top Set images for the real M87* data. We also investigate the 
properties of the image reconstructions when removing individ- 
ual baselines. 


G.1. Synthetic data 


In order to test the efficacy of the imaging pipelines, we generate 
synthetic data sets using both geometric models and snapshots of 
GRMHD simulations. All of the geometric models contain 0.6 Jy 
in compact flux density and 0.5 Jy in a larger extended jet mod- 
eled by three Gaussians. We show the ring model as an example 
in Fig. G.1. The sizes and locations of these three Gaussians are 
identical to those reported in Table 10 of M 87* 2017IV, with the 
flux densities scaled down to total 0.5 Jy rather than 0.6 Jy. Of the 
geometric models we generated, four of them were utilized in the 
procedure to generate the Top Set images. These four geometric 
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models are very similar to those used in M 87* 2017IV and are 

as follows: 

1. cres180: An asymmetric ring model with rọ = 23 uas, a 
brightness position angle oriented south, and blurred by a 
circular Gaussian beam of FWHM 10 pas. 

2. ring: A thin uniform ring of radius rọ = 23 yas blurred by a 
circular Gaussian beam of FWHM 10 pas. 

3. dblsrc: Two circular Gaussian components each with 
FWHM of 20 as. One is located at the origin with a flux 
density of 0.27 Jy, while the second is positioned at Ag A. = 
30 uas and Agel, = — 12 was with a flux density of 0.33 Jy. 

4. disk: A uniform disk of radius rọ = 35 as blurred by a 
circular Gaussian beam of FWHM 10 pas. 

Ground-truth images of these models are shown in Fig. G.2. 
Aside from these training data, we also generated seven valida- 
tion data sets to ensure good performance from the CLEAN and 
RML method Top Sets. These data, along with reconstructions 
from each pipeline, are shown in Fig. G.3. The model descrip- 
tions are as follows 

1. cres-968, cres8, cres90: Three asymmetric ring models 
with ro = 23 uas, and brightness position angles oriented 
east, north, and west, blurred by a circular Gaussian beam of 
FWHM = 10 as. 

2. edisk: A uniform elliptical disk model with a major axis 
of 66 uas, a minor to major axis ratio of 0.65, a major axis 
position angle of 60°, and a blurring of 10 pas. 

3. point+disk: A point source plus symmetric disk model 
containing a 10 as point source centered in 100 as diam- 
eter disk. The point source to disk flux ratio is chosen to be 
0.192. 

4. point+edisk: A point source plus elliptical disk model 
containing a 10 as point source centered in an ellipse of 
major axis of 96 uas, a minor to major axis ratio of 0.8, a 
major axis position angle of 60°, and a blurring of 10 yas. 
The point source to disk flux ratio is 0.16. 

5. GRMHD: A snapshot from a GRMHD simulation of M87*. 
The synthetic — visibilities were generated using 

eht-imaging and all include random thermal noise, station- 

based gain and phase errors, and polarimetric leakage. For a 
more in-depth explanation of synthetic visibility generation, we 
refer readers to Appendix C.2 of M 87* 2017IV. 

Figures G.2 and G.3 show the image reconstructions for 
four training and seven validation data sets including the 
GRMHD model from each imaging pipeline, respectively. The 
image reconstructed with the fiducial parameters is displayed 
for DIFMAP, eht-imaging, and SMILI. A single image ran- 
domly selected from the posterior is displayed for THEMIS and 
Comrade. The THEMIS reconstruction of the GRMHD syn- 
thetic data uses a Hybrid Themage model. This figure clearly 
shows that a very wide range of morphology types can be recov- 
ered by all imaging pipelines. These simulations are based on 
the same (u, v) coverage of the 2018 EHT observation, and thus 
have a common point spread function (PSF). This ensures that 
the recovered structures by imaging pipelines are not artificially 
caused by the (u,v) coverage and/or the PSF, but rather come 
from the real structure imprinted in the data. 

Based on owx values and ring parameters, we again con- 
firmed that most of the Top Set parameters can reconstruct 
the correct morphology. For instance, 95, 96, 97% of cres-98 
April 21 band 3 images reconstructed with the Top Set param- 
eters from DIFMAP, eht-imaging, SMILI passed the pyx crite- 
ria imposed on the training data sets though some images have 
deviated structure as seen in the training data sets. Figure G.4 
shows the significance of the deviated structure in the Top Set of 
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Fig. G.1. Example geometric ring model used to generate the synthetic 
data. The left panel exhibits the morphological characteristics of the 
ring on a linear scale, encompassing a FOV of 130 uas. The right panel 
depicts the logarithmic scale representation of the extended jet model 
(FOV = 2900 uas). 


synthetic data reconstructions with SMILI. Above a certain 
dynamic range threshold (e.g., 10% of peak intensity), artifi- 
cial structures which deviate from the ground truth are shown 
but its fraction in the Top Set is minor, especially the repeated 
central structure. In addition, the fiducial images recover the 
station gains as presented in Fig. G.5. This shows the multi- 
plicative gains at each telescope, in comparison with true gains 
from the cres188 model on April 21 band 3. In a same man- 
ner, the station gains are derived for M87* which are consistent 
across different imaging approaches, as presented in Fig. G.6. 
Certain imaging pipelines may perform poorly at certain epochs 
and bands. In case of point+disk April 21 band 3, only 40% 
of DIFMAP images reconstructed with the Top Set parameters 
passed the pyx criteria. 


G.2. The distribution of M87* Top Set images 


As shown in the main text, all Top Set M87* images have con- 
sistent ring features. However, there are slight differences among 
the Top Set reconstructions. Figure G.7 shows 25 randomly 
selected Top Set reconstructions from HOPS band 3 April 21 
data for DIFMAP, eht-imaging, and SMILI. We see the ~ 40 yas 
ring for all images, and most images have a similar bright- 
ness feature, but the ring width seems to have relatively large 


uncertainty depending on the image assumptions of each imag- 
ing method. Also, a few images from eht-imaging and SMILI 
show double-ring structures although they are a minority in the 
Top Sets. As discussed above, we see similar minor double-ring 
structures in ring and crescent synthetic data reconstructions, 
so the second ring is unlikely to be real and we anticipate this 
does not affect our conclusion. These trends are also seen on 
the other observational dates and bands. Image variations among 
Top Sets are relatively larger, especially for data with poor 
coverage. 


G.3. Image dependence on sites 


In addition to analyzing the dependence of M87* images on fre- 
quency bands and changes in coverage from different epochs, it 
is necessary to understand how much the reconstructed image 
relies on data from single antenna sites. For this purpose, we 
imaged M87* excluding all the data from baselines connected to 
one site, and repeated the imaging for all six sites. The test was 
performed on April 21 band 3 data, using the eht- imaging and 
Comrade pipelines as representative for the RML and Bayesian 
methods. The first row of Fig. G.8 shows the mean image 
reconstructed by Comrade, the second row shows eht-imaging 
images reconstructed with fiducial parameters, while the third 
one presents eht-imaging fiducial images of the 2017 April 
11 data, reconstructed with the 2018 pipeline. We see that in 
most cases, the removal of one antenna in 2018 affects the qual- 
ity of the image reconstruction, both for the Comrade and the 
eht-imaging pipelines. 

The strongest image degradation happens when removing 
the Chile sites, which are fundamental to reconstruct the basic 
image morphology, but also imaging without the SMT or the 
Hawai 5 sites results in the appearance of a tessellation pattern 
for eht- imaging and in a loss of contrast for Comrade. Remov- 
ing PV baselines results in an elongation of the ring-like struc- 
ture for both pipelines, while the absence of either the GLT or the 
LMT stations only introduces minor artifacts, without changing 
the overall image morphology. 

Compared to 2017, the 2018 eht-imaging reconstructions 
without one antenna are affected by a stronger tessellation pat- 
tern. However, the morphology of the ring-like structure is more 
robust to the removal of a single site in 2018 compared to 
2017. For example, in 2018, a ring-like structure is still distin- 
guishable even without the Chile site. Also the removal of the 
LMT antenna does not affect the orientation of the ring, and the 
removal of PV does not strongly affect the morphology. 
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Fig. G.2. Four training geometric models as imaged by each imaging method. The first row shows the visibility amplitudes of the model (orange) 
compared to the visibility amplitudes measured for M87" (blue). The second row shows the ground-truth images. The DIFMAP, eht-imaging, and 
SMILI rows show a fiducial image made from the same parameter sets as the images shown in Fig. 17. The THEMIS and Comrade rows show a 
random draw from the posterior. The circle in the DIFMAP panels represent a Gaussian blurring kernel of 20 was. The white lines in the THEMIS 
and Comrade panels represent the size of FWHM for the Gaussian blurring kernel needed to match the DIFMAP effective resolution. 
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Fig. G.3. Seven geometric validation models plus one GRMHD snapshot as imaged by each method. The first row shows the visibility amplitudes 
of the model (orange) compared to the visibility amplitudes measured for M87* (blue). The second row shows the ground-truth images. The 
DIFMAP, eht-imaging, and SMILI rows show a fiducial image made from the same parameter sets as the images shown in Fig. 17. The THEMIS 
and Comrade rows show a random draw from the posterior. The THEMIS reconstruction of the GRMHD synthetic data uses a Hybrid Themage 
model. The circle in the DIFMAP panels represent a Gaussian blurring kernel of 20 uas. The white lines in the THEMIS and Comrade panels 
represent the size of FWHM for the Gaussian blurring kernel needed to match the DIFMAP effective resolution. 
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Fig. G.4. Stacked images with a dynamic range cutoff used to investigate the structural deviations in the Top Set (from SMILI as a representative). 
First, each Top Set image is normalized with its peak intensity, and the pixels with a flux larger than a certain threshold (e.g., 0.1) are converted to 
unity (unless 0). Then all the images across the Top Set are stacked and normalized by the maximum number of Top Set images so that the fraction 
of reliable features over the entire Top Set is emphasized. The ground-truth structure of each synthetic data model is shown in gray contours. 
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Fig. G.5. Multiplicative gain correction factors at each station as a function of time. The ground truth gains (gray circles) are used to generate the 
synthetic data with the cres180 model for April 21, band 3. Model gains from each pipeline are derived from the fiducial images or posterior 


samples of the cres180 synthetic data reconstructions. 
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Fig. G.6. Multiplicative gain correction factors at each station as a function of time from the fiducual images or posterior samples of the image 


reconstructions for the M87* HOPS band 3 data. 
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Fig. G.7. Twenty five randomly selected Top Set images from HOPS band 3 April 21 data for DIFMAP (upper left), eht- imaging (upper right), 
and SMILI (bottom). 
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Fig. G.8. Example reconstructions of M87* for 2018 April 21 and 2017 April 11 after omitting visibilities to each geographical site. The images 
presented show reconstructions that exclude all baselines to the indicated site (i.e., mimicking an observation without that site). Top and middle 
rows show the reconstructions for Comrade and eht-imaging respectively for band 3 on 2018 April 21. The images for the eht-imaging 
pipeline were reconstructed using the eht - imaging fiducial parameters (see Sect. 5.2). The images shown for Comrade are the mean images from 
the posterior. The bottom row shows the reconstructions for 2017 April 11 using the 2018 eht- imaging pipeline but 2017 fiducial parameters 


(M 87* 2017 IV). The ellipse in each panel shows the corresponding synthesized beam with uniform weighting, but the image is not convolved 
with this beam. 
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Appendix H: Imaging with DIFMAP using a point 
source model for initial self-calibration 


This Appendix reports the results of the DIFMAP parameter sur- 
vey, which used an initial phase-only self-calibration process 
assuming a point source starting model. Unlike the parameter 
survey discussed in the main text (Sect. 5.1.1), which seeks to 
determine the best geometric model for initial self-calibration, 
this survey requires a careful consideration of the CLEAN win- 
dow location. This is due to the fact that the self-calibration pro- 
cess results in the brightest part of the source being placed at the 
map's coordinate origin, requiring a shift in the position of the 
window center's position for the subsequent CLEAN and self- 
calibration procedures. 

To ensure optimal performance, the ideal window center 
location should be specified for each source model and not 
treated as a global free parameter. In this analysis, we employed 
an automated optimal mask search procedure during the param- 
eter survey using a point source model for initial self-calibration. 
The procedure involved a systematic search of the parameter 
space by shifting the window center from the map origin in 
increments of 5 uas along the right ascension and declination 
directions. The shift is limited to a maximum value equal to 
the size of the CLEAN window. To maintain consistency, the 
five imaging parameters presented in Sect. 5.1.1 were kept con- 
stant at their fiducial parameters which were obtained from the 
CLEAN imaging of the 2017 EHT M87* data (M 87* 2017 IV). 
These parameters included a compact flux density of 0.5 Jy, stop- 
ping CLEANing when the required compact flux density was 
reached, an ALMA weight re-scaling factor of 0.1, a window 
diameter of 60 pas, and a (u, v) weight exponent of —1. The pro- 


band 1 band 2 


April 21 


cedure computed the Nee value for each image in the survey, as 
described in Sect. 5.2, and determined the optimal window posi- 
tion based on the image with the lowest x value. This approach 
enabled the determination of the optimal window position for 
each source model and window diameter. 

Next, we conducted a parameter survey by varying the five 
imaging parameters (Sect. 5.1.1), assuming the optimal window 
position obtained during the automatic mask search procedure. 
We followed the image selection criteria outlined in Sect. 5.2 
to select both Top Set and fiducial images from the survey. 
The fiducial images resulting from this survey are displayed in 
Fig. H.1. Our findings indicate that the survey was able to accu- 
rately reconstruct the majority of the synthetic data models. For 
the real M87* data, the Top Set images displayed a ring-like mor- 
phology, consistent with the results from the survey discussed 
in Sect. 5.1.1. This result demonstrates that the ring-like struc- 
ture for M87*, reconstructed by DIFMAP, is stable independently 
of the geometrical model selected for the first phase of self- 
calibration. However, the overall performance of this survey was 
lower than that of the survey that optimized the geometric mod- 
els for initial phase-only self-calibration, resulting in a smaller 
number of Top Set images. This can be attributed to the fact that 
the point source model used as initial phase-only self-calibration 
model in this survey is overly simplistic for many geometric 
models and actual M87* source structure. Despite encounter- 
ing challenges with the uniform disk model, resulting in a lower 
number of Top Set images obtained from the survey, success- 
ful reconstruction of both the ring and crescent synthetic models 
was achieved. These models were found to closely resemble the 
true source structure. 
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Fig. H.1. Fiducial images produced by the DIFMAP pipeline reconstructed through the implementation of an initial phase-only self-calibration 
process using a point source model. The reconstructed images exhibit ring-like structures, which are comparable to the results obtained from the 
pipeline that utilized closure phases to select the optimal geometric model for the initial phase-only self-calibration, as demonstrated in Fig. 7. The 
circle represents a Gaussian blurring kernel with a FWHM of 20 jas. 
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Appendix I: Testing the backward compatibility of M87* data, including the real data from the HOPS pipeline 
the DIFMAP pipelines on the 2017 data (April 11, low-band) and the synthetic data sets for the train- 
Gu . ing. Our analysis revealed that the resulting images align well 
The DIFMAP pipelines presented in Sect. 5.1.1 are updated as — with those previously reported in (M 87* 2017 IV). Our fiducial 
compared to the DIFMAP pipelines used in previous EHT imag- images of the M87* data and synthetic data can be found in 
ing analyses (M 87* 2017IV; Sgr A* 2017 III). We conducted a Fig. I.1. 
backward compatibility test of these pipelines on the 2017 EHT 
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Fig. I.1. Fiducual images of the 2017 April 11 low-band data obtained from the 2018 DIFMAP pipelines discussed in Sect. 5.1.1, as a backwards 
compatibility check. The top row of images is obtained using the pipeline that conducts initial phase-only self-calibration by employing a point 
source model, whereas the bottom row of images is obtained from the pipeline that identifies the optimal geometric model for self-calibration. The 
circles represent a Gaussian blurring kernel with a FWHM of 20 jas. 
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Appendix J: THEMIS imaging model details and 
priors 


The THEMIS image reconstructions are composed of three dif- 
ferent primary model components, the splined raster model, a 
large-scale asymmetric Gaussian, and a thin slashed crescent. 
A more detailed description of the large-scale Gaussian and 
crescent models, as well as details about their technical imple- 
mentation, are described in Broderick et al. (2020a). A much 
more detailed discussion of the raster model can be found in 
Broderick et al. (2020b). 

The splined raster model is defined by a rectilinear set of 
control points which may independently vary in intensity, 7j y. 
The intensity map is produced using an approximate cubic spline 
interpolation between the control points. The field of view FOV, 
and FOV, are permitted to vary, along with the orientation of 
the grid $ and the raster center (x, y). The large-scale Gaussian 
is asymmetric, and is characterized by its total flux Jg, a sym- 
metric standard deviation o'g, asymmetry parameter Ag, position 
angle $c, and position (xc, yo). The slashed crescent is based on 
the xs-ringauss model also described in Broderick et al. (2020b), 
which also serves as the core of the THEMIS xs-ringauss geo- 
metric model featured in Sect. 6.2 and in M 87* 2017 VI. 

As a part of the Hybrid Themage model, we eliminate the 
contributions from the asymmetry and pinned Gaussian parame- 
ters by appropriately restricting their priors. We also restrict the 
priors on the fractional width parameter y to range from 0% to 
596. This ensures the ring width is well below that of the raster 
spacing, to avoid excessive flux trading between the ring and 
raster components. This leaves the ring flux Jy, outer radius Rout, 
linear brightness gradient fx, brightness gradient position angle 
ox, and position (xx, yx) as searchable parameters. 

We use a deterministic even-odd (DEO) swap tempering 
scheme with a Hamiltonian Monte-Carlo sampling kernel based 
on the Stan package. Chain convergence is assessed based on 
traditional criteria such as the split R, auto-correlation time, and 
visual inspection of the parameter traces. The parameters and 
priors for each model component are listed in Table J.1. 


Table J.1. THEMIS Raster Model Priors. 


Component Parameter Units Prior? 
Raster ImN Jyuas? — £(10714,3 x 103) 
FOV, pas U0, 200) 
FOV, pas U0, 200) 
rad U(-0.25z, 0.257) 
x pas U(—40, 40) 
y pas U(—40, 40) 
Gaussian Ig Jy U0, 10) 
OG mas U(0.1, 10^) 
Ag U0, 1) 
óc rad U(0, 7) 
XG mas U(-2, 2) 
ya mas U(-2, 2) 
Ring Ix Jy U(0, 2) 
Rout pas U0, 107) 
y i35 U(0, 0.05) 
fx " U0, 1) 
oy rad U(-7, n) 
Xx pas U(—40, 40) 
yx pas U(—40, 40) 


è Logarithmic priors from a to b are represented by L(a, b), and uniform 
(linear) priors from a to b are represented by U(a, b). 


A79, page 57 of 63 


The Event Horizon Telescope Collaboration: A&A, 681, A79 (2024) 


Appendix K: THEMIS Raster dimensions 


Table K.1. THEMIS raster size (N, x Ny) survey for April 21 Band 3. 


Raster Dims. red. y? Alog(Z) 
4x4 1.15 -6 
5x5 1.14 - 
6x6 1.17 -21 
7x7 1.20 -20 


The best fit raster dimensions N, x N, for the THEMIS Raster 
models in principle depend on the (u,v) coverage and source 
complexity. To construct the best model, one should survey the 
raster dimensions for each unique data set. In order to provide 
easier comparisons across data sets, we only survey the raster 
dimensions on our primary science data set, April 21 band 3, 
and use those best fit raster dimensions for the other data sets. 
We present the results of this survey in Table K.1, where we sur- 
vey a 4 X 4, 5 x 5, 6 x 6, and 7 x 7 raster. We rank the models in 
terms of the difference in the logarithm of the Bayesian Evidence 
(Alog(Z)), where positive values are more preferred. The 5 x 5 
raster is most preferred by this data, followed by the 4 x 4 raster. 
The 5 x 5 raster is also the same dimensions as the model used 
in Broderick et al. (20222), which applied the Hybrid Themage 
model to the 2017 data. Thus, the image structures we show in 
this analysis should be directly comparable to the image struc- 
tures in the previous paper. We also report the best fit reduced 
X? for the surveyed raster sizes, and note that all models provide 
good fits to the data. 


Table K.2. Complex visibility reduced x? (x7) for the THEMIS Raster 
and Hybrid Themage model best fits to the HOPS pipeline data. 


THEMIS Raster Hybrid Themage 
Red. y; Red. yy, Alog(Z) 
April 21 band 1 1.41 1.43 -9 
band 2 1.1 1.17 +1 
band 3 1.14 1.07 +26 
band 4 1.19 1.14 +23 
April 25 band 1 2.20 2.7 +6 
band 2 1.71 2.05 +8 
band 3 1.17 1.17 +18 
band 4 0.96 0.97 +14 


We can also compare the fit quality between the standard 
THEMIS raster and Hybrid Themage models by similarly com- 
paring the Alog(Z). In Table K.2 we show the reduced y? and 
difference in Bayesian evidence between the Hybrid Themage 
and raster-only model for both days and all bands. Here we see 
that the Hybrid Themage model is generally preferred in nearly 
all data sets, except for the April 21 band 1 data. 


Appendix L: Comrade image model and prior 


For all Comrade image reconstructions, we construct our image 
model in the following steps. We begin with a rasterized image 
model convolved with a B-spline kernel of order 3 (see Eqs. L.1 
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and L.2) to generate flux densities that smoothly vary in all direc- 
tions. We set the raster points by (/;,mj;) given by Eq. L.3, for a 
given FOV and number of pixels N, and N, for each side. The 
FOV and number of pixels are the hyperparameters of the model: 


l-l m-mj 
I(l,m) = ) B| — |B lij, L.1 
(l,m) D(a B p (L.1) 
where B(x) is a B-spline kernel of order 3 given by, 
x?/2 O<sx<l 
B(x) = $ (—2x? + 6x -—3)/2 1<x<2 (L.2) 
(3 — x)*/2 2<x<3 
and raster points are given by, 
li = -FOV /2 + s,/2 + s,(i- 1 
[2 + Sx/2 + sx(i— 1) (L3) 


mj = —FOV/2 + Sy/2 ks syCj us 1) 


where s, — FOV/N, and s, — FOV/N,. We restrict the images 
considered in this paper to having N, = Ny, such that the 
images are square and s, — s,. We add a circular Gaussian with 


FWHM -2 V21nec: = 1 mas (Gimas) to the raster to model the 
amplitudes and station gains of short intra-site (ALMA-APEX, 
JCMT-SMA) baselines. We form our final image model, 


M(ij, f, fe) = fü zs te) x Kj) + fof X Gimas 


where f is the total flux density of the raster (also described 
as compact flux) and f, is the fraction of the total flux density 
that corresponds to the flux density of the Gaussian. We take 
the Fourier transform of M(I;, f, fe) to get model Fourier trans- 
form 7, for each baseline ab. The Fourier amplitudes are then 
corrupted to model the individual station gain amplitudes |g| to 
construct the visibility amplitudes, Iul = lgallgollZaol 

We form the visibility amplitude likelihood and closure 
phase likelihood as described in Appendix F of Sgr A* 2017 IV. 
Next, we form our prior distribution by using a uniform distri- 
bution U(0.0, 1.5) Jy for the prior on the total flux density f and 
U(0.0, 1.0) for the fraction of the total flux density f, for the flux 
density of the large-scale Gaussian. For the raster pixels fluxes, 
Ij, we use a symmetric Dirichlet distribution 


(KS) T ea 
pile) = Ét, 
rex l | 


(L.4) 


(L.5) 


where I is the flattened vector of pixel fluxes, K is the total 
number of pixels in the image, I is gamma distribution, and 
é is the concentration parameter that controls the sparsity and 
smoothness of the image. The Dirichlet support is on the sim- 
plex, meaning that >); I; = 1 and 0 x J; < 1. For this work, we 
set £ — 1, which is equivalent to a uniform distribution on the 
simplex. 

From the analysis of the crossing baseline tracks and the 
EHT station flux density calibration parameters (Table A.1), 
for station gain priors, we used a Normal distribution for the 
log-gain amplitudes for each station. For the April 21 data, we 
imposed prior widths of 10% or all stations, that is a N(0.0, 0.1) 
prior, except 30% and 100% prior widths for LMT and GLT, 
respectively. We used the same gain priors for the April 25 data 
except for PV which used a 100% prior width, since it was found 
that PV exhibited large amplitude fluctuations after UTC 02:00 
due to poor weather as mentioned in Appendix D. 
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Appendix M: Impact of compact flux density on 
image reconstruction 


The non-imaging constraints on compact flux density permit a 
wide range of allowable values, and a comparison of the com- 
pact flux densities measured by the Bayesian methods suggests 
that the current data set cannot be used to constrain the compact 
flux on horizon scales. In this appendix, we test the robustness 
of the ring reconstructions from the RML and CLEAN methods 
against different assumptions of the compact flux density. 
Firstly, we investigated the dependence of the assumed total 
compact flux density (Fp) on the estimated ring parameters 
of M87* from reconstructed images following the approach 
adopted in our previous EHT analysis (M 87* 2017 IV). In this 
analysis, we only change Fepct, while all other imaging parame- 
ters are kept at their fiducial values. We select acceptable images 
based on two criteria: (1) normalized Xp» aie ca € 2 (with a 
10% systematic uncertainty for DIFMAP), and (2) a correspond- 
ing self-calibration solution with 0.9 < median(1/|gswrl) < 1.1, 
where gsmr is the gain for the SMT station. Figure M.1 illus- 
trates how the estimated ring parameters vary with the assumed 
total compact flux density. While the ring width for all pipelines 
and fractional central brightness for RML methods slightly 
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increase with Fepct, the ring diameter and position angle val- 
ues of the acceptable images do not exhibit significant depen- 
dence on F;y4. Therefore, the ring diameter and position angle 
are robust against different assumptions on the compact total flux 
density. 

Secondly, We tested if our DIFMAP imaging script could dis- 
tinguish different morphologies for different values of the total 
compact flux density using synthetic data. We prepared a series 
of synthetic data with the same morphologies of four geomet- 
ric models used for Top Set selection in the imaging survey, but 
with the assumed compact total flux density in range from 0.4 
to 1.1Jy. Here, we used the total flux density, which is the sum 
of the assumed compact total flux density and flux density of 
a larger extended jet modeled by three Gaussians, to be 1.1Jy 
for all synthetic data. We conduct image reconstruction using 
the DIFMAP imaging script with the fiducial parameter sets, but 
adopting the compact flux densities to be the same as that of 
the assumed values of synthetic data. In Fig. M.2, we show the 
reconstructed images of geometric models with different values 
of the total compact flux density. It is clearly demonstrated that 
our DIFMAP imaging script can distinguish ring-like morpholo- 
gies from disk or double source morphologies if we properly 
choose the expected total compact flux density. 
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Fig. M.1. Measured diameter d, width w, position angle, and fractional central brightness fc of M87*, measured from image reconstructions assum- 
ing different total compact flux densities, Fp. for DIFMAP (orange), eht-imaging (blue), and SMILI (green) (see Sect. 5.3). All measurements 
are made using REx. All other imaging parameters were set to the fiducial parameters of the corresponding pipeline. DIFMAP values were measured 
after restoring with a 20 uas FWHM Gaussian beam. The solid lines indicate the measured value, and the shaded regions give the +1o uncertainty 
of the REx measurement. The darker colored regions correspond to values of F. that produces images that have (1) normalized Xip Xog ca <2 
(with 10% systematic uncertainty for DIFMAP), and (2) a corresponding self-calibration solution with 0.9 < median(1/|gsmr|) < 1.1. 
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dblsrc ring cres180 


disk 


Fig. M.2. Reconstructed DIFMAP images of the four geometric models with different values of the total compact flux density. We use the fiducial 
parameters from the 2018 Top Set, but with different assumed total compact flux densities. The top label above each column indicates the assumed 
total compact flux density in units of Jy. 
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Appendix N: Geometric modeling details 


* xs — ringauss 
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0 1/2 ln 3/20 2m 

Fig. N.1. Schematic summary of the xs-ringauss and mF-ring model 
parameters. The xs-ringauss is defined to be a disk with a circular hole 
removed, that is allowed to be offset from the center of the disk and 
rotated. A circular Guassian is pinned to the center of the disk as an 
emission floor, and an additional elliptical Gaussian is pinned to the 
inner edge of the offset hole. The mF-ring is an m-ring of order 2 that 
has an elliptical stretch which is allowed to be rotated. A flat emission 
floor is included to match the interior of the m-ring. Both models have 
an additional blur defined on them with a Gaussian blurring kernel. The 
intensity profiles of the fiducial models along two axes are indicated in 
red and green. An angular intensity profile is also shown for the mF-ring 
in blue. 


Image reconstructions of M87* show strong support for a ring- 
like structure. We quantify this support through the Bayesian evi- 
dence as a robust metric for the goodness of fit. 


The Bayesian approach is reliant on making an appropri- 
ate estimate of the posterior distribution, P(O|M, D), for some 
parameters O of a chosen model, M, that is conditioned on some 
data, D. The posterior is defined by the relationship, 


P(D|M, ©) P(9|M) " L£(O|M) x(0|M) 


FOI M.D) = 777 PDM ZM) ^ 


(N.D 


which also defines the likelihood of the data for a chosen model 
and parameter, 


£(O|M) = PDM, 9) (N.2) 
the prior probability on model parameters, 

m(@|M) = P(O|M) (N.3) 
and the Bayesian evidence for a chosen model, 

Z(M) = P(D|M) = [ cow 1(0|M) dO. (N.4) 


The Bayesian evidence allows the relative support between two 
models to be determined through a Bayes factor, 

Z(M1) 7M2) 

Z(M1) 7M2) 


R(™M,, Mb) = (N.5) 


where z(M) is the probability prior defined on the set of models. 
We use the Bayesian evidence with a constant prior to define the 
goodness of fit for a model M € {M} as, 

A log(Z) = log Z(M) - maximum [log Z( {M} )], (N.6) 
where the argument of Alog(Z) is the Bayes factor between M 
and the best performing model in the set. 

The best performing models are selected with this metric 
to construct fiducial models for feature comparison; namely, a 
stretched m-ring of order 2 with a flat emission floor and two 
elliptical Gaussians (mF-ring), and an xs-ringauss with two ellip- 
tical Gaussians an a large scale circular Gaussian. The schematic 
summary of defining parameters of these models is shown in 
Fig. N.1. See Sect. 6.2 for details on these models. 

The shared geometric parameters of interest between both 
the mF-ring and the xs-ringauss model are their diameters, 
widths, central brightnesses, and brightness asymmetry position 
angles. The diameter of the mF-ring is defined to be the sum of 
the debiased semi-major and semi-minor axis, while that of the 
xs-ringauss is taken to be the sum of its inner and outer radii, 


qc fo + Îi, 
Rout + Rin, 


Here, the debiased radius (7;) is linked to the defining radius of 
the mF-ring (7;), and the FWHM of the blurring Gaussian kernel 
(c7) through: 


mF-ring (N.D) 
xs-ringauss. 


g*2 


f;-rj- = N.8 
aT Fé InQ)r; Oe) 
We define the fractional widths of the mF-ring and the xs- 
ringauss respectively to be, 

Te 


Rou—Rinto* 


mF-ring 
fw = (N.9) 


, Xs-ringauss. 
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Table N.1. THEMIS xs-ringauss model parameters and priors. 


Parameter Description Units Prior? 
Io+h Flux Density of Crescent and Fixed Gaussian Jy U(0, 2) 
Rout Outer Radius of Crescent pas U0, 100) 
dy Position Angle rad U(-7, r) 
Fractional Thickness U(0.0001, 0.9999) 
T Structural Asymmetry U(0.0001, 0.9999) 
fx Flux Asymmetry U(0, 1) 
g* Width of Gaussian Smoothing Kernel Has *4(0, 100) 
Ip Flux Density of Central Gaussian Emission Floor Jy U(0, 2) 
Of Width of Central Gaussian Emission Floor pas U(0, 25) 
Ig Flux Density of Large-Scale Gaussian Jy U0, 10) 
OG Width of Large-Scale Gaussian mas L0 (72, 10) 
I, /Ip Ratio of Fixed Gaussian Flux Density to Crescent Flux Density U0, 1) 
Ox/Rout Ratio of Fixed Gaussian Width to Outer Radius of Crescent U0, 3) 
Oy/o~x Axis Ratio of Fixed Gaussian *4(0, 100) 
I, Flux Density of Additional Gaussian Component Jy «4(0, 2) 
Xo Central x-coordinate of Additional Gaussian pas U(-200, 200) 
yo Central y-coordinate of Additional Gaussian Has U(—200, 200) 
Tg Width of Additional Gaussian pas U0, 100) 
Ag Anisotropy of Additional Gaussian U0, 0.99) 
[n Position Angle of Additional Gaussian rad U(0, 7) 


a Logarithmic priors in base x spanning a range x^ to x^ are represented by £,(a, b), and uniform (linear) priors spanning a to b are represented by 


Ua, b). 


Table N.2. mF-ring parameters and priors. 


Component Parameter Description Units Prior 

mF-ring ri Semi-major axis length Has U(10, 30) 
ro/T| Semi-major/Semi-major length ratio *4(0.1, 1) 
Pm Semi-major axis orientation rad U (0,7) 
fe Relative central brightness U (0.0, 1.0) 
Ai Amplitude of the i® mode U (0.0, 0.5) 
Bi Phase of the it" mode U (0.0, 0.5) 
o*/(2-V2In@)) Blurring kernel standard deviation uas — 4(1.0,50) 
f mF-ring/(all nuisance Gaussians) flux ratio U (0, 1) 

Nuisance Je Primary/secondary Gaussian flux ratio T U(0.51, 1) 

Gaussians C gi Standard deviation of i" nuisance Gaussian Has U(2, 40) 
Ti Semi-major/semi-minor length ratio of i* nuisance Gaussian vu *4(0.1, 1) 
gi Stretch angle of i™ nuisance Gaussian rad U(0, 7) 
Xi Horizontal displacement of primary Gaussian from m-ring center pas 14(—50, 50) 
x2 Horizontal displacement of secondary Gaussian from primary center pas 14(—50, 50) 
MI Vertical displacement of primary Gaussian from m-ring center Has 14(—50, 50) 
y2 Vertical displacement of secondary Gaussian from primary center Has 14(—50, 50) 

and the brightness asymmetry position angle to be, which is given by, 
Eu brightness emission foor , mF-ring 
à = b +m, m Tng N.10) fc _ em 7 Ko n (N.11) 
$x, xs-ringauss. brightness of central Gaussian xs-ringauss. 


The last shared parameter shown between the two models is the 
fractional central flux depression, or relative central brightness, 
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As described above in the main text, the THEMIS xs- 
ringauss model used to fit the data from the 2018 observations 
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Fig. N.2. Representative images for each observing band from random posterior samples for the Hybrid Themage, xs-ringauss, and mF-ring 


models, for the April 21 and April 25 HOPS data. 


is the same model used to fit the data from the 2017 observa- 
tions. We list the model parameters and their priors in Tables N.1 
and N.2. 

Even though we construct and fit the Hybrid Themage, xs- 
ringauss, and mF-ring in the visibility domain, we can also pro- 
duce conventional image domain representations of the mod- 
els. Figure N.2 shows the on-sky representations of the Hybrid 
Themage, xs-ringauss, and mF-ring models. We find that these 
models reproduce many of the same features seen in the more 
agnostic models shown in Sect. 5, such as a consistent diam- 


eter and position angle of the brightest part of the ring in 
the southwest. The xs-ringauss and mF-ring models use a pair 
of asymmetric Gaussians as nuisance components, designed 
to help fit the residual data that cannot be captured by the 
ring component. The raster component in the Hybrid Themage 
serves a similar purpose. We note that for April 21 band 3 
and band 4, the non-ring components settle in the west and 
southwest part of the image, but these components occupy 
more varied locations in the band 1, band 2, and April 25 
reconstructions. 
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